René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"Incorrect and Inconsistent output of function PSD when scale_by_freq=False","articleBody":"If we state `scale_by_freq=False`, then we are asking for the power spectrum of the signal, instead of the power spectral density; i.e. units of V^2 instead of V^2/Hz, for example. \n\n`mlab.psd`, which actually is `mlab._spectral_helper` under the hood, does not return the expected value of the power in a sinusoid. Furthermore, in contrast to the PSD case, it does not correct for the energy of the signal lost to the window function, meaning that the results differ depending on the window used. \n\nHere's an illustration. If we have a simple signal of a cosine oscillation happening at 100Hz with amplitude of 10, we would expect the power spectrum result at 100Hz to be 50, regardless of the window used. (A cosine like this that is identically zero at each integer second means we can use one second FFTs and compare the rectangular window to the hanning window, without worrying about any distortion)\n\n In MATLAB:\n\n```\nfs=2048;\ntt=[0:1/fs:60-1/fs];\nx = 10*cos(2*pi*100*tt);\n[pR, f]=pwelch(x, ones(1,fs), 0, fs, fs, 'power');\n[pH, f]=pwelch(x, hann(fs), 0, fs, fs, 'power');\n```\n\ngives `pR(101)==50` and `pH(101)==50.000...`, as expected. \n\nIn `matplotlib.mlab`:\n\n```\nfs = 2048\ntt = np.arange(fs*64)/float(fs)\nx = 10*np.cos(2*np.pi*100*tt)\n# One second FFT segments\npsdHann, ff = mlab.psd(x, Fs=fs, NFFT=fs, window=mlab.window_hanning, scale_by_freq=False)\npsdRect, ff = mlab.psd(x, Fs=fs, NFFT=fs, window=mlab.window_none, scale_by_freq=False)\n```\n\ngives `psdRect[100] == 102400`, which too large by a factor of `nfft` and `psdHann[100] ==  68233.3`, which differs further by some window correction factor. \n\nI came across this in scipy/scipy#4682, where I was adapting the code for inclusion in scipy's signal processing repertoire. There, I used `scale = 1.0 / win.sum()**2` for the power spectrum case. I haven't looked for the exact places where `_spectral_helper` should be fixed, at this time. \n","author":{"url":"https://github.com/e-q","@type":"Person","name":"e-q"},"datePublished":"2015-04-12T02:18:33.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":3},"url":"https://github.com/4328/matplotlib/issues/4328"}