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?