pythonscipysignal-processingfftspectrogram

Scipy Incorrect Amplitude when computing FFT


When trying to compute the FFT of this stacked Sinus Signal the resulting FFT does NOT match the initial input. The code initializes a result array that is populated with sinus signals of varying frequency, but constant amplitude. Later I try to re-create what the STFT function does (which btw. also yields the same result). Now the problem is that choosing a 2^n window size results in periodical errors in the FFT Amplitude. Some even more than 15%. Why is that so? All info I could find online was to specifically use a 2^n window size. When choosing an arbitrary window size e.g. 10000 the problem does not occur...

import numpy as np
from scipy.fft import fft, fftshift

# Define the sample rate and the duration of the signal
sample_rate = 10000  # in Hz, common sample rate in audio
duration = 100     # in seconds, 10 milliseconds

# Create a time array
t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)

# Frequencies of the sine waves
frequencies = np.arange(1000,4000, 100)
result = np.zeros(t.shape[0])

# Create sine waves
for freq in frequencies:
    result += np.sin(2 * np.pi * freq * t)

def compute_fft_sections(data, sample_rate, x):
    n = len(data)
    section_length = 8192
    overlap = int(section_length * (0.5))

    fft_results = []
    normalised_sections=[]
    for i in range(0, n - overlap, overlap):

        # Stop if the end of the array is reached
        if i + section_length >= n:
            break

        section = data[i:i + section_length]

        # Apply Hanning window
        window = np.hanning(len(section))
        windowed_section = section * window

        # Normalize the windowed section
        # Here we normalize by the sum of the window values
        normalized_section = windowed_section / window.sum()
        normalised_sections.append(section)
        normalised_sections.append(normalized_section)

        # fft + downsample symetric
        fft_result = fft(normalized_section)
        fft_result = fftshift(fft_result)
        fft_result = fft_result[:len(fft_result)//2]
        #fft_result = fft_result[::len(fft_result)]

        # append to results
        fft_results.append(fft_result*2)

    return fft_results, normalised_sections

# Example usage
fft_sections, normalised_sections = compute_fft_sections(result,
                                    sample_rate=0.0001, 
                                    x=120)

I tried varying the window width to achieve correct amplitudes. But for some reason using non-2^n wide windows reduces the error... which is exactly the opposite I expected.

Periodical SciPy FFT Error:

Periodical SciPy FFT Error


Solution

  • Why is that so? All info I could find online was to specifically use a 2^n window size.

    The reason you see this advice isn't for the sake of accuracy - it's for the sake of speed.

    The FFT works by recursively splitting the window, and recombining the answers. For that reason, you get the best performance using a number which is a power of 2. You get average performance if you pick a number which is a product of factors 2, 3, and 5, such as any power of 10. You get the worst performance if you pick a prime window size.

    FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, by using symmetries in the calculated terms. The symmetry is highest when n is a power of 2, and the transform is therefore most efficient for these sizes. For poorly factorizable sizes, scipy.fft uses Bluestein’s algorithm and so is never worse than O(n log n).

    (Source.)

    SciPy’s FFTPACK has efficient functions for radix {2, 3, 4, 5}, so this returns the next composite of the prime factors 2, 3, and 5 which is greater than or equal to target. (These are also known as 5-smooth numbers, regular numbers, or Hamming numbers.)

    (Source.)

    When choosing an arbitrary window size e.g. 10000 the problem does not occur...

    The window size of 10000 has a nice property - it contains frequencies that are the exact frequencies you're using to produce the sine waves.

    In contrast, if you look at the frequencies obtained from a window size of 8192, all of them are slightly different from the frequencies used to generate the signal.

    Example:

    import matplotlib.pyplot as plt
    
    freq = scipy.fft.fftfreq(8192, 1 / 10000)
    deviation = []
    for frequency in frequencies:
        closest = freq[np.abs((freq - frequency)).argmin()]
        deviation.append(closest - frequency)
        print(frequency, "closest freq", closest)
    plt.scatter(frequencies, deviation)
    plt.xlim([0, 5000])
    

    This produces a plot where the X axis represents the frequency of a sine wave in Hertz, and the Y axis represents the difference between the sine wave frequency and the nearest FFT frequency. From this, you can see that the two differ by as much as 0.6 Hz.

    plot of frequency versus frequency deviation

    This is a plot of the FFT, using the same window length.

    fft

    Notice that when the frequencies match, the FFT magnitude is big, but when the frequencies don't match, the FFT magnitude is small.