soundtouch change rate matlab implementation

Posted 2021-01-21 fellow1988

soundtouch implement of changing rate in a way same with resample(SRC).

 

%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);