Porting code from MATLAB

Question

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

Answer 1

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}]