diff --git a/src/audio/kaiser_window.pl b/src/audio/kaiser_window.pl deleted file mode 100755 index 1a2e50bfa..000000000 --- a/src/audio/kaiser_window.pl +++ /dev/null @@ -1,210 +0,0 @@ -#!/usr/bin/perl -w - -use warnings; -use strict; - -# The resampling algorithm: https://ccrma.stanford.edu/~jos/resample/ -# https://www.mathworks.com/help/signal/ref/kaiser.html -# "Thus kaiser(L,beta) is equivalent to -# besseli(0,beta*sqrt(1-(((0:L-1)-(L-1)/2)/((L-1)/2)).^2))/besseli(0,beta)." -# Matlab kaiser calls besseli(): -# https://www.mathworks.com/help/matlab/ref/besseli.htm -# https://en.wikipedia.org/wiki/Bessel_function - -sub print_table { - my $tableref = shift; - my $name = shift; - my @table = @{$tableref}; - my $comma = ''; - my $count = 0; - print("static const float $name = {\n "); - foreach (@table) { - print("$comma$_"); - #print(sprintf("%.6f\n", $_)); - if (++$count > 4) { - $count = 0; - print(",\n "); - $comma = ''; - } else { - $comma = ', '; - } - } - print("\n};\n\n"); -} - - -use POSIX (); - -# This is a "modified" bessel function, so you can't use POSIX j0() -sub bessel { - my $x = shift; - - my $i0 = 1; - my $f = 1; - my $i = 1; - - while (1) { - my $diff = POSIX::pow($x / 2.0, $i * 2) / POSIX::pow($f, 2); - last if ($diff < 1.0e-21); - $i0 += $diff; - $i++; - $f *= $i; - } - - return $i0; -} - -sub kaiser { - my $L = shift; - my $beta = shift; - my @retval; - - #print("L=$L, beta=$beta\n"); exit(0); - - for (my $i = 0; $i < $L; $i++) { - my $val = bessel($beta * sqrt(1.0 - - POSIX::pow( - ( - ( - ($i-($L-1.0)) - ) / 2.0 - ) / (($L-1)/2.0), 2.0 )) - ) / bessel($beta); - - unshift @retval, $val; - } - return @retval; -} - - -my $zero_crossings = 5; -my $bits_per_sample = 16; -my $samples_per_zero_crossing = 1 << (($bits_per_sample / 2) + 1); -my $kaiser_window_table_size = ($samples_per_zero_crossing * $zero_crossings) + 1; - -# if dB > 50: 0.1102 * ($db - 8.7) -my $db = 80.0; -my $beta = 0.1102 * ($db - 8.7); - -my @table = kaiser($kaiser_window_table_size, $beta); - -print_table(\@table, 'kaiser_window'); - -# Kaiser window has "sinc function" ("cardinal sine") applied to it: -# sin(pi * x) / (pi * x) -# "For example, to use the ideal lowpass filter, the table would contain -# h(l) = sinc(l/L)." - -use Math::Trig ':pi'; -for (my $i = 1; $i < $kaiser_window_table_size; $i++) { - my $x = $i / $samples_per_zero_crossing; - $table[$i] *= sin($x * pi) / ($x * pi); -} - -print_table(\@table, 'with_sinc'); - -# "Our implementation also stores a table of differences ¯h(l) = h(l + 1) − h(l) between successive -# FIR sample values in order to speed up the linear interpolation. The length of each table is -# Nh = LNz + 1, including the endpoint definition ¯h(Nh) = 0." - -my @differences = (); -for (my $i = 1; $i < $kaiser_window_table_size; $i++) { - push @differences, $table[$i] - $table[$i - 1]; -} -push @differences, 0; - -print_table(\@differences, 'differences'); - - -# Might as well use this code as a test harness... - -use autodie; -my $fnamein = shift @ARGV; -my $fnameout = shift @ARGV; -my $inrate = shift @ARGV; -my $outrate = shift @ARGV; - -print("Resampling $fnamein (freq=$inrate) to $fnameout (freq=$outrate).\n"); - -open(IN, '<:raw', $fnamein); -my @src = (); - -# this assumes mono Sint16 raw data since we aren't parsing .wav files. -# !!! FIXME: deal with multichannel audio. -my $channels = 1; - -# this is inefficient, but this is just throwaway code... -while (read(IN, my $bytes, 2) == 2) { - my ($samp) = unpack('s', $bytes); - push @src, $samp; -} - -close(IN); - -my $ratio = $outrate / $inrate; -my $sample_frames_in = scalar(@src) / $channels; -my $sample_frames_out = $sample_frames_in * $ratio; - -my $outsamples = $sample_frames_out * $channels; -#my @dst = (0) x ($outsamples); -my @dst = (); -print("Resampling $sample_frames_in input frames to $sample_frames_out output (ratio=$ratio).\n"); - - -my $inv_spzc = int(POSIX::ceil(($samples_per_zero_crossing * $inrate) / $outrate)); -my $padding_len; -if ($ratio < 1.0) { - $padding_len = int(POSIX::ceil(($samples_per_zero_crossing * $inrate) / $outrate)); -} else { - $padding_len = $samples_per_zero_crossing; -} - -# You need to pad the input or we'll get buffer overflows. -# !!! FIXME: deal with multichannel audio. -for (my $i = 0; $i < $padding_len; $i++) { - push @src, 0; - unshift @src, 0; -} - -# !!! FIXME: deal with multichannel audio. -my $time = 0.0; -for (my $i = 0; $i < $outsamples; $i++) { - my $srcindex = int($time * $inrate); # !!! FIXME: truncate or round? - - my $ftime = $srcindex / $inrate; # this would be $time if we didn't convert $srcindex to int. - my $fnexttime = ($srcindex + 1) / $inrate; - - # do this twice to calculate the sample, once for the "left wing" and then same for the right. - my $sample = 0; - my $interpolation = 1.0 - ($fnexttime - $time) / ($fnexttime - $ftime); - my $filterindex = int($interpolation * $samples_per_zero_crossing); - - $srcindex += $padding_len; - - for (my $j = 0; ($filterindex + ($j * $samples_per_zero_crossing)) < $kaiser_window_table_size; $j++) { - $sample += int($src[$srcindex - $j] * ($table[$filterindex + $j * $samples_per_zero_crossing] + $interpolation * $differences[$filterindex + $j * $samples_per_zero_crossing])); - } - - $interpolation = 1 - $interpolation; - $filterindex = $interpolation * $samples_per_zero_crossing; - for (my $j = 0; ($filterindex + ($j * $samples_per_zero_crossing)) < $kaiser_window_table_size; $j++) { - $sample += int($src[$srcindex + 1 + $j] * ($table[$filterindex + $j * $samples_per_zero_crossing] + $interpolation * $differences[$filterindex + $j * $samples_per_zero_crossing])); - } - - push @dst, $sample; - - # "After each output sample is computed, the time register is incremented by 2nl+nη /ρ (i.e., time is incremented by 1/ρ in fixed-point format)." - $time += 1.0 / $outrate; -} - -open(OUT, '>:raw', $fnameout); - -# this is inefficient, but this is just throwaway code... -foreach (@dst) { - print OUT pack('s', $_); -} - -close(OUT); - -print("Done.\n"); -