#! /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";