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