soundtouch implement of changing rate in a way same with resample(SRC).
When rate < 1, it need interpolate sample. and delete samples when rate > 1.
After interpolation, there may be introduce high frequency. To avoid aliase, we usally apply a low pass filter on interpolated signal.
For the case of deleting samples, some signal may contains high frequecy, the difference between sample may be very sharp. Therefore, we use low pass filter before deleting samples.
When design fir low pass filter, we usaully use "w" (rad/s) as the main parameter instead freqence "f".
The relationship of "w" and "f" as following:
w = 2*pi * f/ fs;
For cutoff frequency, wc = 2*pi * fc / fs = 2 * pi * fc/(2*fn) = pi * fc / fn, where fc is cutoff frequency, fn is Nywuist frequency.
For fc = = fn, wc = pi;
We use cubic method to interplate sample, the principle of cubic interpolation refer to https://www.paulinternet.nl/?page=bicubic
%calc low pass filter coefficient. The low pass filter based on sinc function with hamming window.
function coeff = calCoeffs(cutoffFreq, len)
coeff = zeros(len ,1);
wc = 2 * pi * cutoffFreq;
tempCoeff = 2 * pi / len;
sum = 0;
for i = 0 : 1 : len -1
cntTemp = (i - len/2);
temp = cntTemp * wc;
% sinc function
if temp ~=0
h = sin(temp) / temp;
else
h = 1;
end
%hamming window
w = 0.54 + 0.46 * cos(tempCoeff * cntTemp);
coeff(i+1) = w * h;
sum = sum + coeff(i+1);
end
coeff = coeff / sum;
end
function output = firfilter(input, coeff)
inputLen = length(input(:, 1));
filterLen = length(coeff(:, 1));
output = zeros(inputLen ,1 );
outputLen = inputLen - filterLen;
for i = 1: 1: outputLen
inpos = i;
sum = 0;
for j = 1:1:filterLen
sum = sum + input(inpos ,1) * coeff(j, 1);
inpos = inpos + 1;
end
output(i, 1) = sum;
end
end
function output = cubicInterpolation(input, rate)
inputLen = length(input(:,1));
outputLen = floor(inputLen / rate);
output = zeros(outputLen ,1);
inputIdx = 1;
fract = 0;
outputIdx = 1;
while inputIdx < inputLen - 4
x1 = fract;
x2 = x1 * x1;
x3 = x1 * x2;
p0 = input(inputIdx , 1);
p1 = input(inputIdx + 1 , 1);
p2 = input(inputIdx + 2, 1);
p3 = input(inputIdx + 3, 1);
output(outputIdx ,1) = (-0.5*p0 + 1.5*p1 -1.5 *p2 + 0.5*p3) * x3 +(p0 - 2.5*p1 + 2*p2 -0.5*p3) *x2 + (-0.5*p0 + 0.5*p2) * x1 + p1;
outputIdx = outputIdx + 1;
fract = fract + rate;
whole = floor(fract);
fract = fract - whole;
inputIdx = inputIdx + whole;
end
end
function output = linearInterpolation(input, rate)
inputLen = length(input(:,1));
outputLen = floor(inputLen / rate);
output = zeros(outputLen ,1);
inputIdx = 1;
fract = 0;
outputIdx = 1;
while inputIdx < inputLen - 4
p0 = input(inputIdx , 1);
p1 = input(inputIdx + 1 , 1);
output(outputIdx ,1) = (1-fract) * po + fract * p1;
outputIdx = outputIdx + 1;
fract = fract + rate;
whole = floor(fract);
fract = fract - whole;
inputIdx = inputIdx + whole;
end
end
function output = changeRate(input, rate, interpMethod)
inputLen = length(input(:, 1));
outputLen = floor(inputLen / rate);
output = zeros(outputLen, 1);
if rate > 1
cutoffFreq = 0.5 / rate;
else
cutoffFreq = 0.5 * rate;
end
filterLen = 64;
coeff = calCoeffs(cutoffFreq, filterLen);
if rate < 1
%slow down, need interpolation first;
if strcmp(interMethod, 'cubic')
output = cubicInterpolation(input, rate);
else
output = linearInterpolation(input, rate);
end
output = firfilter(output, coeff);
else
%fast, need filter out the high freqency, then delete samples
output = firfilter(input, coeff);
if strcmp(interMethod, 'cubic')
output = cubicInterpolation(output, rate);
else
output = linearInterpolation(output, rate);
end
end
end
main.m:
clc;
clear all;
[input fs] = wavread('input.wav');
%if do SRC, rate = inputfs / outputfs;
rate = 0.5;
output = changeRate(input, rate, 'cubic');
wavwrite(output, fs, 'output.wav);