| #! /usr/bin/perl -w |
| # |
| # DirectSound |
| # |
| # Copyright 2011-2012 Alexander E. Patrakov |
| # |
| # This library is free software; you can redistribute it and/or |
| # modify it under the terms of the GNU Lesser General Public |
| # License as published by the Free Software Foundation; either |
| # version 2.1 of the License, or (at your option) any later version. |
| # |
| # This library is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| # Lesser General Public License for more details. |
| # |
| # You should have received a copy of the GNU Lesser General Public |
| # License along with this library; if not, write to the Free Software |
| # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA |
| |
| use strict; |
| use Math::Trig; |
| |
| # This program generates an array of Finite Impulse Response (FIR) filter |
| # values for use in resampling audio. |
| # |
| # Values are based on the resampler from Windows XP at the default (best) |
| # quality, reverse engineered by saving kvm output to a wav file. |
| |
| # Controls how sharp the transition between passband and stopband is. |
| # The transition bandwidth is approximately (1 / exp_width) of the |
| # Nyquist frequency. |
| |
| my $exp_width = 41.0; |
| |
| # Controls the stopband attenuation. It is related but not proportional |
| # to exp(-(PI * lobes_per_wing / exp_width) ^2) / lobes_per_wing |
| |
| my $lobes_per_wing = 28; |
| |
| # Controls the position of the transition band and thus attenuation at the |
| # Nyquist frequency and above. Amended below so that the length of the FIR is |
| # an integer. Essentially, this controls the trade-off between good rejection |
| # of unrepresentable frequencies (those above half of the lower of the sample |
| # rates) and not rejecting the wanted ones. Windows XP errs on the side of |
| # letting artifacts through, which somewhat makes sense if they are above |
| # 20 kHz anyway, or in the case of upsampling, where we can assume that the |
| # problematic frequencies are not in the input. This, however, doesn't match |
| # what linux resamplers do - so set this to 0.85 to match them. 0.98 would |
| # give Windows XP behaviour. |
| |
| my $approx_bandwidth = 0.85; |
| |
| # The amended value will be stored here |
| my $bandwidth; |
| |
| # The number of points per time unit equal to one period of the original |
| # Nyquist frequency. The more points, the less interpolation error is. |
| my $fir_step = 120; |
| |
| |
| # Here x is measured in half-periods of the lower sample rate |
| sub fir_val($) |
| { |
| my ($x) = @_; |
| $x *= pi * $bandwidth; |
| my $s = $x / $exp_width; |
| my $sinc = $x ? (sin($x) / $x) : 1.0; |
| my $gauss = exp(-($s * $s)); |
| return $sinc * $gauss; |
| } |
| |
| # Linear interpolation |
| sub mlinear($$$) |
| { |
| my ($y1, $y2, $mu) = @_; |
| return $y1 * (1.0 - $mu) + $y2 * $mu; |
| } |
| |
| # to_db, for printing decibel values |
| sub to_db($) { |
| my ($x) = @_; |
| return 20.0 * log(abs($x))/log(10.0); |
| } |
| |
| my $wing_len = int($lobes_per_wing / $approx_bandwidth * $fir_step + 1); |
| $bandwidth = 1.0 * $lobes_per_wing / $wing_len; |
| |
| my $amended_bandwidth = $bandwidth * $fir_step; |
| my $fir_len = 2 * $wing_len + 1; |
| my @fir; |
| |
| # Constructing the FIR is easy |
| for (my $i = 0; $i < $fir_len; $i++) { |
| push @fir, fir_val($i - $wing_len); |
| } |
| |
| # Now we have to test it and print some statistics to stderr. |
| # Test 0: FIR size |
| print STDERR "size: $fir_len\n"; |
| |
| # Test 1: Interpolation noise. It should be less than -90 dB. |
| |
| # If you suspect that 0.5 is special due to some symmetry and thus yields |
| # an abnormally low noise figure, change it. But really, it isn't special. |
| my $testpoint = 0.5; |
| |
| my $exact_val = fir_val($testpoint); |
| my $lin_approx_val = mlinear($fir[$wing_len], $fir[$wing_len + 1], |
| $testpoint); |
| |
| my $lin_error_db = to_db($exact_val - $lin_approx_val); |
| |
| printf STDERR "interpolation noise: %1.2f dB\n", $lin_error_db; |
| |
| # Test 2: Passband and stopband. |
| # The filter gain, ideally, should be 0.00 dB below the Nyquist |
| # frequency and -inf dB above it. But it is impossible. So |
| # let's settle for -80 dB above 1.08 * f_Nyquist. |
| |
| my $sum = 0.0; |
| $sum += $_ for @fir; |
| |
| # Frequencies in this list are expressed as fractions |
| # of the Nyquist frequency. |
| my @testfreqs = (0.5, 0.8, 1.0, 1.08, 1.18, 1.33, 1.38); |
| foreach my $testfreq(@testfreqs) { |
| my $dct_coeff = 0.0; |
| for (my $i = 0; $i < $fir_len; $i++) { |
| my $x = 1.0 * ($i - $wing_len) / $fir_step; |
| $dct_coeff += $fir[$i] * cos($x * $testfreq * pi); |
| } |
| printf STDERR "DCT: %1.2f -> %1.2f dB\n", |
| $testfreq, to_db($dct_coeff / $sum); |
| } |
| |
| # Now actually print the FIR to a C header file |
| |
| chdir ".." if -f "./make_fir"; |
| open FILE, ">dlls/dsound/fir.h"; |
| select FILE; |
| |
| print "/* generated by tools/make_fir; DO NOT EDIT! */\n"; |
| print "static const int fir_len = $fir_len;\n"; |
| print "static const int fir_step = $fir_step;\n"; |
| print "static const float fir[] = {\n"; |
| |
| for (my $i = 0; $i < $fir_len; $i++) { |
| printf "%10.10f", $amended_bandwidth * $fir[$i]; |
| if ($i == $fir_len - 1) { |
| print "\n"; |
| } elsif (($i + 1) % 5 == 0) { |
| print ",\n"; |
| } else { |
| print ", "; |
| } |
| } |
| print "};\n"; |