How can I vectorize the calculation of a multi-taper spectrum?

Question

Normally I compute the spectrum of a signal using pmtm:

signal = rand(1000,1);
NW = 4;
Fr = 1:50;
Fs = 200;
[p, fr] = pmtm( signal, NW, Fr, Fs);

However I'm looking for a way to vectorize this so I can compute multiple spectra at the same time. I tried:

signal = rand(1000,10); %<--- notice I have 10 columns instead of 1
NW = 4;
Fr = 1:50;
Fs = 200;
[p, fr] = pmtm( signal, NW, Fr, Fs);

but it produces an error that doesn't really tell me what I did wrong. I know I can wrap the call to pmtm in a loop.

Here is the error:

Error using .* Matrix dimensions must agree.

Error in pmtm>mtm_spectrum (line 231) [Xx,w] = computeDFT(E(:,1:k).*x(:,ones(1,k)),nfft,Fs);

Error in pmtm (line 142) [S,k,w] = mtm_spectrum(x,params);

This leads me to suspect that there isn't a vectorized way to achieve what I want. I was hoping that someone here would know how to do this.

`mtm_spectrum` forces `x` to be vector in the line `x = x(:);`. If `x` is matrix, it will have times n_col larger dimension and this is why you get matrix dimension error. Have you tried arrayfun? It's pseudo-vectorized solution, but may speed up the execution time a little. Something like `arrayfun(@(x)pmtm(signal(:,x), NW, Fr, Fs),1:size(signal,2))`. But you can get only one output variable with this. — yuk, Mar 03 '13 at 23:27
You probably can't avoid a for loop here, which is quite expensive because `pmtm` is slow. If performance is your concern, then you can speed it up significantly by precalculating the slepian sequences with `dpss` and pass them on to `pmtm`. — Max, Sep 17 '16 at 13:30

Lolo · Answer 1 · 2013-03-04T05:35:50.227

What makes you think that pmtm is designed to work on matrices? The help menu specifically expects a vector:

>> help pmtm
 pmtm   Power Spectral Density (PSD) estimate via the Thomson multitaper 
    method (MTM).
    Pxx = pmtm(X) returns the PSD of a discrete-time signal vector X in 
    the vector Pxx.  Pxx is the distribution of power per unit frequency.
    The frequency is expressed in units of radians/sample.  pmtm uses a 
    default FFT length equal to the greater of 256 and the next power of
    2 greater than the length of X.  The FFT length determines the length
    of Pxx.
...

If you are after performance, you may want to get the spectral representation of your signal with a 2D fft (fft2) and then compute the distribution of power on your own. This will allow you to work on 2D arrays without any for loop but that will mean more coding for you for sure :-(

How can I vectorize the calculation of a multi-taper spectrum?

1 Answers1