Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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;

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; 


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

Any suggestions or pointers to improve it?

EDIT: Example usage


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

Elapsed time is 0.305124 seconds.


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

share|improve this question
Can you give an example of input that is suitable for benchmarking? – C. E. Jul 23 '14 at 1:03
@Pickett I added a couple of examples – Wolfy Jul 23 '14 at 1:37
up vote 8 down vote accepted

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

        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}]
share|improve this answer
+1. The difference is probably in your approximation of HammingWindow. – RunnyKine Jul 23 '14 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 '14 at 2:49

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.