Potential Discrepancy in DOLfYN `fft_frequency` for Half Frequencies

[ssolson]

@jmcvey3 can you let me know if I am approaching this incorrectly?

Description:

Issue with Half Frequencies in fft_frequency Function I think there is potential issue in the fft_frequency function, specifically related to how the half frequencies (positive frequencies) are returned. The current implementation might be including an extra frequency bin which should not be present.

https://github.com/MHKiT-Software/MHKiT-Python/blob/2ccc286a65685e5e2f5ab68a467209917f0161d9/mhkit/dolfyn/tools/fft.py#L6-L31

Current Behavior:

The function, when called with full=False, returns np.abs(f[1:int(nfft / 2. + 1)]). This slice appears to include one frequency bin more than expected for the positive frequency range.

Expected Behavior:

I believe the correct implementation for returning half frequencies should be np.abs(f[1:int(nfft / 2.0)]), excluding the additional bin at the end.

Test Case Issue:

In my test case for this function, I encountered a size mismatch error when comparing the returned half frequencies to the expected positive frequencies. This issue arises in the line assert_allclose(freq_half, positive_freqs), where freq_half is expected to match the size and values of positive_freqs.

If the current returned value is correct could you let me know if I am slicing the positive and negative values from the full frequencies incorrectly?

def test_fft_frequency(self):
    fs = 1000  # Sampling frequency
    nfft = 512  # Number of samples in a window

    # Test for full frequency range
    freq_full = tools.fft.fft_frequency(nfft, fs, full=True)
    assert_equal(len(freq_full), nfft)

    # Check symmetry of positive and negative frequencies, ignoring the zero frequency
    positive_freqs = freq_full[1:int(nfft / 2)]
    negative_freqs = freq_full[int(nfft / 2) + 1:]
    assert_allclose(positive_freqs, -negative_freqs[::-1])

    # Test for half frequency range
    freq_half = tools.fft.fft_frequency(nfft, fs, full=False)
    assert_equal(len(freq_half), int(nfft / 2) - 1)
    assert_allclose(freq_half, positive_freqs)  # Ignore the zero frequency

[jmcvey3]

positive_freqs = freq_full[1:int(nfft / 2)] should be positive_freqs = freq_full[1:int(nfft / 2)+1], no? Otherwise you'll be cropping out the Nyquist frequency.

[ssolson]

freq_full[1:int(nfft / 2)+1]

In the example I give if I define "positive_freqs" as suggested my last frequency is negative:

        fs = 1000  # Sampling frequency
        nfft = 512  # Number of samples in a window
        freq_full = tools.fft.fft_frequency(nfft, fs, full=True)
        positive_freqs = freq_full[1:int(nfft / 2)+1]
        print(positive_freqs )

I get

array([   1.953125,    3.90625 ,    5.859375,    7.8125  ,    9.765625,
         ... ,
        490.234375,  492.1875  ,  494.140625,  496.09375 ,  498.046875,
       -500.      ])

Should the last frequency in the "positive frequencies" be negative?

or should positive frequencies be defined as: positive_freqs =f[:int(nfft / 2) + 1] ?

[jmcvey3]

No. I'm not sure why the Nyquist frequency here (500) is negative, but I assume that's why np.abs() is used in the original function

[jmcvey3]

It just appears to be how the numpy function works. The middle frequency technically should be equivalent whether + or -, and we throw out 0 Hz because it's irrelevant for us.

>>> np.fft.fftfreq(10)
array([ 0. ,  0.1,  0.2,  0.3,  0.4, -0.5, -0.4, -0.3, -0.2, -0.1])

If you're concerned about the frequency vector of the calculated PSDs, we can write up a test comparing this code to scipy's Welch function. I compared them at one time and found they gave the same output.