René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"[Bug]: Possible problem with parameter `NFFT` and `pad_to` of `mlab._spectral_helper`","articleBody":"### Bug summary\r\n\r\nIf my understanding is correct:\r\n- `NFFT` of `mlab._spectral_helper` corresponds to `nperseg` of spectral funtions of `scipy.signal`.\r\n  It's the actual length of signal (segment) for DFT.\r\n- `pad_to` of `mlab._spectral_helper` corresponds to `nfft` of spectral funtions of `scipy.signal`.\r\n  It's the length of DFT, i.e. parameter  `n` of `fft` funtion.\r\n  https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L385\r\n\r\nBut the implementation of Matplotlib is really confusion:\r\n- If length of signal is less than NFFT, signal will pad to length of NFFT.\r\n  https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L342-L351\r\n  \u003c!--https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L316-L318--\u003e\r\n  And the window is determined by NFFT instead of actual length of signal.  This would cause wrong window correction.\r\n  https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L376-L380\r\n\r\n- The Nyquist frequency is `NFFT/2` instead of `pad_to/2`.\r\n  https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L411-L416\r\n\r\n### Code for reproduction\r\n\r\n```python\r\nimport numpy as np\r\nfrom matplotlib import mlab\r\nfs = 1000\r\nt = np.arange(0, 1, 1/fs)\r\nrng = np.random.default_rng(42)\r\ns = np.sin(2 * np.pi * 100 * t) + 0.5 * rng.normal(size=t.shape)\r\n\r\nnfft = 2048\r\nnsignal = len(s) #1000\r\npsd, freqs = mlab.psd(s, NFFT=nfft, Fs=fs, window=mlab.window_none)\r\n\r\nmag = np.abs(np.fft.rfft(s, n=nfft))\r\npsd1 = (mag)**2/nsignal/fs\r\nif not nfft % 2:\r\n  psd1[1:-1] *= 2\r\nelse: \r\n  psd1[1:] *= 2\r\nrelative_error = np.abs( 2 * (psd-psd1)/(psd + psd1) )\r\nprint(relative_error.max())\r\n```\r\n\r\n### Actual outcome\r\n\r\n0.6876640419947717\r\n\r\n### Expected outcome\r\n\r\n0","author":{"url":"https://github.com/gapplef","@type":"Person","name":"gapplef"},"datePublished":"2022-12-27T05:46:50.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":5},"url":"https://github.com/24822/matplotlib/issues/24822"}