tools/make_fir

   1 #! /usr/bin/perl -w
   2 #
   3 # DirectSound
   4 #
   5 # Copyright 2011-2012 Alexander E. Patrakov
   6 #
   7 # This library is free software; you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation; either
  10 # version 2.1 of the License, or (at your option) any later version.
  11 #
  12 # This library is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 # Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with this library; if not, write to the Free Software
  19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
  20
  21 use strict;
  22 use Math::Trig;
  23
  24 # This program generates an array of Finite Impulse Response (FIR) filter
  25 # values for use in resampling audio.
  26 #
  27 # Values are based on the resampler from Windows XP at the default (best)
  28 # quality, reverse engineered by saving kvm output to a wav file.
  29
  30 # Controls how sharp the transition between passband and stopband is.
  31 # The transition bandwidth is approximately (1 / exp_width) of the
  32 # Nyquist frequency.
  33
  34 my $exp_width = 41.0;
  35
  36 # Controls the stopband attenuation. It is related but not proportional
  37 # to exp(-(PI * lobes_per_wing / exp_width) ^2) / lobes_per_wing
  38
  39 my $lobes_per_wing = 28;
  40
  41 # Controls the position of the transition band and thus attenuation at the
  42 # Nyquist frequency and above. Amended below so that the length of the FIR is
  43 # an integer. Essentially, this controls the trade-off between good rejection
  44 # of unrepresentable frequencies (those above half of the lower of the sample
  45 # rates) and not rejecting the wanted ones. Windows XP errs on the side of
  46 # letting artifacts through, which somewhat makes sense if they are above
  47 # 20 kHz anyway, or in the case of upsampling, where we can assume that the
  48 # problematic frequencies are not in the input. This, however, doesn't match
  49 # what linux resamplers do - so set this to 0.85 to match them. 0.98 would
  50 # give Windows XP behaviour.
  51
  52 my $approx_bandwidth = 0.85;
  53
  54 # The amended value will be stored here
  55 my $bandwidth;
  56
  57 # The number of points per time unit equal to one period of the original
  58 # Nyquist frequency. The more points, the less interpolation error is.
  59 my $fir_step = 120;
  60
  61
  62 # Here x is measured in half-periods of the lower sample rate
  63 sub fir_val($)
  64 {
  65     my ($x) = @_;
  66     $x *= pi * $bandwidth;
  67     my $s = $x / $exp_width;
  68     my $sinc = $x ? (sin($x) / $x) : 1.0;
  69     my $gauss = exp(-($s * $s));
  70     return $sinc * $gauss;
  71 }
  72
  73 # Linear interpolation
  74 sub mlinear($$$)
  75 {
  76     my ($y1, $y2, $mu) = @_;
  77     return $y1 * (1.0 - $mu) + $y2 * $mu;
  78 }
  79
  80 # to_db, for printing decibel values
  81 sub to_db($) {
  82     my ($x) = @_;
  83     return 20.0 * log(abs($x))/log(10.0);
  84 }
  85
  86 my $wing_len = int($lobes_per_wing / $approx_bandwidth * $fir_step + 1);
  87 $bandwidth = 1.0 * $lobes_per_wing / $wing_len;
  88
  89 my $amended_bandwidth = $bandwidth * $fir_step;
  90 my $fir_len = 2 * $wing_len + 1;
  91 my @fir;
  92
  93 # Constructing the FIR is easy
  94 for (my $i = 0; $i < $fir_len; $i++) {
  95     push @fir, fir_val($i - $wing_len);
  96 }
  97
  98 # Now we have to test it and print some statistics to stderr.
  99 # Test 0: FIR size
 100 print STDERR "size: $fir_len\n";
 101
 102 # Test 1: Interpolation noise. It should be less than -90 dB.
 103
 104 # If you suspect that 0.5 is special due to some symmetry and thus yields
 105 # an abnormally low noise figure, change it. But really, it isn't special.
 106 my $testpoint = 0.5;
 107
 108 my $exact_val = fir_val($testpoint);
 109 my $lin_approx_val = mlinear($fir[$wing_len], $fir[$wing_len + 1],
 110         $testpoint);
 111
 112 my $lin_error_db = to_db($exact_val - $lin_approx_val);
 113
 114 printf STDERR "interpolation noise: %1.2f dB\n", $lin_error_db;
 115
 116 # Test 2: Passband and stopband.
 117 # The filter gain, ideally, should be 0.00 dB below the Nyquist
 118 # frequency and -inf dB above it. But it is impossible. So
 119 # let's settle for -80 dB above 1.08 * f_Nyquist.
 120
 121 my $sum = 0.0;
 122 $sum += $_ for @fir;
 123
 124 # Frequencies in this list are expressed as fractions
 125 # of the Nyquist frequency.
 126 my @testfreqs = (0.5, 0.8, 1.0, 1.08, 1.18, 1.33, 1.38);
 127 foreach my $testfreq(@testfreqs) {
 128     my $dct_coeff = 0.0;
 129     for (my $i = 0; $i < $fir_len; $i++) {
 130         my $x = 1.0 * ($i - $wing_len) / $fir_step;
 131         $dct_coeff += $fir[$i] * cos($x * $testfreq * pi);
 132     }
 133     printf STDERR "DCT: %1.2f -> %1.2f dB\n",
 134         $testfreq, to_db($dct_coeff / $sum);
 135 }
 136
 137 # Now actually print the FIR to a C header file
 138
 139 chdir ".." if -f "./make_fir";
 140 open FILE, ">dlls/dsound/fir.h";
 141 select FILE;
 142
 143 print "/* generated by tools/make_fir; DO NOT EDIT! */\n";
 144 print "static const int fir_len = $fir_len;\n";
 145 print "static const int fir_step = $fir_step;\n";
 146 print "static const float fir[] = {\n";
 147
 148 for (my $i = 0; $i < $fir_len; $i++) {
 149     printf "%10.10f", $amended_bandwidth * $fir[$i];
 150     if ($i == $fir_len - 1) {
 151         print "\n";
 152     } elsif (($i + 1) % 5 == 0) {
 153         print ",\n";
 154     } else {
 155         print ", ";
 156     }
 157 }
 158 print "};\n";