2012-05-06 18:06:13 +00:00
|
|
|
#! /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
|
|
|
|
|
2020-03-09 00:30:22 +00:00
|
|
|
open FILE, ">", "fir.h";
|
2012-05-06 18:06:13 +00:00
|
|
|
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";
|