make_fir 4.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
#! /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";