|
@@ -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");
|
|
|
-
|