5
$\begingroup$

I am doing some massive number crunching in MATLAB which involves millions of PSD estimations. Each data segment has length 41. So I have been using the multi-taper method with nfft=41 (details below in the code). My last simulation took more than eight hours. The code is optimized and vectorized as I can make it so I thought maybe if I change nfft=64 that might speed up. The profiler confirms that pmtm is the slowest part by far.

So, looking at the MATLAB help pages, I see that in the documentation it says nothing more than if nfft is larger than the length of the data segment then it is zeropadded. My question is what exactly is MATLAB doing? How exactly is the signal padded?

I can't really find anything googling and I ran my own experiments and I am not getting what I'm supposed to. Here I present my code and the result.


close all
clear all

dseg = [ 395.9640
  392.3630
  388.8230
  385.2820
  381.6130
  378.2600
  374.7470
  371.1860
  367.9920
  364.7830
  361.6630
  358.5660
  355.5620
  352.5820
  349.5530
  346.8110
  343.8130
  340.9010
  338.1460
  335.5110
  332.6300
  330.0470
  327.4020
  324.0930
  321.4680
  318.8760
  316.4470
  314.0850
  311.5310
  309.0600
  306.7570
  304.6690
  302.5820
  300.3890
  298.2480
  296.3650
  294.1660
  291.7610
  290.3300
  287.9700
  285.9850];

% Original signal with length=41
[ffty1,fftf1] = pmtm(dseg,4,41,1/30);
semilogy(fftf1,ffty1,'-o')
hold on

% Matlab "zeropads" to length=64
[ffty2,fftf2] = pmtm(dseg,4,64,1/30);
semilogy(fftf2,ffty2,'-or')

% Original zeropaded on both sides
dsegpad = [zeros(11,1); dseg; zeros(11,1); 0];
[ffty3,fftf3] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf3,ffty3,'-ok')

% Original zeropadded at the end
dsegpad = [dseg; zeros(23,1)];
[ffty4,fftf4] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf4,ffty4,'-og')

% Original signal just repeated
dsegpad = [dseg; dseg(1:23)];
[ffty5,fftf5] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf5,ffty5,'-oc')

I compare the original PSD with MATLAB's padded PSD and they are close. Then I zeropad myself instead of letting MATLAB do it. I zeropadded on both sides of the signal. I tried zeropadding all on one side. I even tried periodically repeating the signal to pad the length. But in all three of the cases my padded PSDs are way larger than the original and MATLABs PSDs. Here is what I see.

MATLAB PSD comparison

I tried looking at pmtm code but got lost pretty fast. Anybody knows what MATLAB is doing and why does it obfuscate it so?

|improve this question|||||
$\endgroup$
  • 1
    $\begingroup$ If you look inside pmtm you will see that the input is transformed by some discrete prolate spheroidal sequences before the zeropadding takes place. That may explain the weird curves you see. I have no idea how you can make the pmtm call run faster. Maybe you can benchmark that function alone to get some idea if its worth writing your own speed-optimized pmtm function. $\endgroup$ – niaren May 2 '13 at 10:41
  • $\begingroup$ @niaren: Not sure what's going on here, but the Kaiser window is supposed to be an approximation of the DPSS window that isn't as computationally expensive. ccrma.stanford.edu/~jos/sasp/Kaiser_Window.html ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html Also I'm just stabbing in the dark here, but is it spending time calculating the exact same DPSS window every time it's called, and that could be pulled out of the function and made into a single call? $\endgroup$ – endolith May 2 '13 at 19:56
  • $\begingroup$ @FixedPoint, anything missing on my answer? $\endgroup$ – Royi Aug 31 '19 at 18:51
1
$\begingroup$

Lets say you have a vector $ x = {\left[ 1, 2, 3, 4 \right]}^{T} $.
You want to have a look on its DFT transform then you apply DFT on it and have the 4 points DFT transform of the data.

In MATLAB it would be something like:

vXDft = fft(vX);

Yet, what happens if you want the 64 point DFT of the same data?
In MATLAB it would be something like:

vXDft = fft(vX, 64);

So what happens?
MATLAB just adds zero elements at the end of the vector until its size is 64:

vXDft = fft([vX; ones([60, 1])]);

This will match what happens above (I assumed Column Vector).

In practice you get a denser grid in the Frequency Domain (64 Points).
Where the points are interpolated by Dirichlet Kernel (The Discrete equivalent of the Sinc function).

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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