# Porting code from MATLAB

I'd like to switch from MATLAB to Mathematica for my research in digital signal processing. I used Mathematica a lot during the '90s but I haven't used it since then, so I might as well be a complete newbie. As a first experiment I tried to port the following MATLAB function (a slow implementation of CQT):

function cq= slowQ(x, minFreq, maxFreq, bins, fs) % x must be a row vector

Q = 1/(2^(1/bins)-1);
maxK = ceil(bins*log2(maxFreq/minFreq));
maxN = ceil(Q*fs/minFreq);
if size(x,2) < maxN
x(end+1:maxN) = 0;
end

for k=1:maxK
fk = minFreq * 2^((k-1)/bins);
N = ceil(Q*fs/fk);
cq(k) = x(1:N) * (hamming(N) .* exp( -2*pi*i*Q*(0:N-1)'/N)) / N;
end

end


This is the code I came up with, which not only looks horrible but is also horribly slow (several orders of magnitude slower than MATLAB):

slowCQT[x_, minFreq_, maxFreq_, bins_, fs_] :=
Module[{Q, maxK, maxN, k, fk, Ncq, cqt},
Q = 1/(2^(1/bins) - 1);
maxK = Ceiling[bins*Log2[maxFreq/minFreq]];
maxN = Ceiling[Q*fs/minFreq];
If[Length[x] < maxN,
x = Join[x, ConstantArray[0, maxN - Length[x]]]];
cqt = ConstantArray[0, maxK];
For[k = 1, k < maxK, k++,
fk = minFreq*2^((k - 1)/bins);
Ncq = Ceiling[Q*fs/fk];
cqt[[k]] =
N[x[[1 ;; Ncq]].(N[HammingWindow[Range[-1/2, 1/2, 1/(Ncq - 1)]]] *
N[Exp[-2*Pi*I*Q*Range[0, Ncq - 1]/Ncq]])/Ncq];];
cqt]


Any suggestions or pointers to improve it?

EDIT: Example usage

MATLAB

fs = 44100;
T = (0:2*fs)/fs;
x=sin(2*pi*440*T);
tic; cqt = slowQ(x,27.5,880,24,fs); toc

Elapsed time is 0.305124 seconds.


Mathematica

fs = 44100;
T = N[Range[0, 2*fs]/fs];
x = N[Sin[2 Pi 440 * T]];
tic = AbsoluteTime[];
cqt = slowCQT[x, 27.5, 880, 24, fs];
AbsoluteTime[] - tic

89.755638

-
Can you give an example of input that is suitable for benchmarking? –  Pickett Jul 23 at 1:03
@Pickett I added a couple of examples –  Wolfy Jul 23 at 1:37

Bear in mind that I'm not completely sure what you are calculating, here is my 5 cents towards what can help you get a better performance. I haven't benched marked the MATLAB code on my system, but the changes implemented in Mathematica lead to a runtime decrease from 78.036 s down to to 0.171 s.

The slowing factors where mainly that you handle a lot of the mathematics with arbitrary precision math (which is slow) only to then afterwards convert to numerical results. If you convert to numerics initially then the operations are carried out much much faster. So rather than doing N@f[Range[0,100]] you can do f[Range[0,100,1.], where 1. is simply such that the range is numerical precision. In addition to this, I changed out the window function since it seemed rather slow, this might be a case where the built in function is over engineered compared to your needs, but I'm not certain of this. You can verify that they are similar using Plot[{hamW[x], HammingWindow[x]}, {x, -1, 1}].

In addition to this I switched out the For loop with a Table. I'm not sure if this impacts performance but I prefer the readability.

So here's the code:

hamW[v_] = 21./46 Cos[2.Pi v] + 25./46

newCQT[x_, minFreq_, maxFreq_, bins_, fs_] :=
Module[{Q, maxK, maxN, k, fk,Ncq, cqt,xn},

Q = 1/(2^(1/bins) - 1);
maxK = Ceiling[bins*Log2[maxFreq/minFreq]];
maxN = Ceiling[Q*fs/minFreq];
xn = N@If[Length[x] >= maxN,x,Join[x, ConstantArray[0, maxN - Length[x]]]];
cqt = ConstantArray[0, maxK];

Table[
fk = minFreq*2^((k - 1)/bins);
Ncq = Ceiling[Q*fs/fk];
xn[[1 ;; Ncq]].(hamW[Range[-0.5, 0.5, 1./(Ncq - 1)]] *
N[Exp[-2*Pi*I*Q*Range[0, Ncq - 1,1.]/Ncq]])/Ncq,{k,1,maxK}]
]


Update: The initial code had an issue due to a poorly chosen window function approximation which is solved now thanks to RunnyKine.

The benchmark was carried out using:

fs = 44100;
T = N[Range[0, 2*fs]/fs];
x = N[Sin[2 Pi 440 * T]];
First@AbsoluteTiming[oldResult = slowCQT[x, 27.5, 880, 24, fs]]
(* 78.036029 *)

First@AbsoluteTiming[newResult = newCQT[x, 27.5, 880, 24, fs]]
(* 0.171415 *)

Max@Abs[newResult - oldResult]
(* 1.25761*10^-15 *)

ListLinePlot[{Log@Abs@oldResult, Log@Abs@newResult}]

-
+1. The difference is probably in your approximation of HammingWindow. –  RunnyKine Jul 23 at 2:45
@RunnyKine Thank you, yes, a minor change to the window function solved the issue perfectly. With the largest error now being on the order of 10^-15. –  jVincent Jul 23 at 2:49