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.
I tried looking at pmtm
code but got lost pretty fast. Anybody knows what MATLAB is doing and why does it obfuscate it so?