3
$\begingroup$

I don't understand why FFT return different maximum amplitude as the signal length increase. I would except that with a large signal length according to the frequency, detected amplitude will be very accurate.

Plot amplitude from FFT against signal length

From where this periodic signal is coming ?

Here is the related python code I used to generate the plot:

%matplotlib qt
%load_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt

def get_fft(x, dt):
    n = len(x)
    fft_output = np.fft.rfft(x)
    rfreqs = np.fft.rfftfreq(n, d=dt)
    fft_mag = [np.sqrt(i.real ** 2 + i.imag ** 2) / n for i in fft_output]
    return np.array(fft_mag), np.array(rfreqs)

def build_signal(amp, freq, signal_len):

    f = freq
    A = amp
    dt = 1
    t = signal_len

    x = np.arange(0, t, dt)
    y = A * np.cos(2*np.pi*f*x)

    fft_mag, rfreqs = get_fft(y, dt)
    return x, y, fft_mag, rfreqs, amp, freq, signal_len

# Build sin wave from 1 to 5000 signal length with freq=1e-3Hz and amplitude=0.5
all_t = []
all_amp = []
for t in np.arange(1, 5000, 100):
    x, y, fft_mag, rfreqs, amp, freq, signal_len = build_signal(amp=0.5, freq=0.001, signal_len=t)

    all_t.append(t)
    all_amp.append(fft_mag.max())

# Plot amplitude from FFT against signal length
plt.figure()
plt.plot(all_t, all_amp, 'o-')
plt.xlabel("Signal length")
plt.ylabel("Amplitude from FFT (0.25 is excepted)")

As MackTuesday ask in comment, I plotted the signal and related FFT for signal length = 1000 and 1400.

signal length = 1000 signal length = 1400

$\endgroup$
  • $\begingroup$ You're right that it shouldn't look this way. I don't see a problem in your code, assuming rfft() does the right thing. It would be interesting to see a comparison of the entire magnitude spectra for signal_len=1000 and signal_len=1400 or so. $\endgroup$ – MackTuesday May 22 '14 at 17:14
  • $\begingroup$ So maximum amplitude accuracy does not seem to be related to signal length but more related to how well our signal is periodic for the specific period. $\endgroup$ – HadiM May 22 '14 at 20:15
5
$\begingroup$

Your code is bit unclear, especially generation of your signal. Python allows for vectorized operations so it is good to use it. What's more, it is good to clearly specify the sampling frequency of your signal and use it then. Also please remember to normalize your FFT by length of your signal (in this particular case) and multiply by 2 (half of spectrum is removed so energy must be preserved).

Phenomena you are facing is connected with fact that your frequency is not matching exactly frequency bin for it. Because of that energy is leaking into other frequency bins. This is due to fact that you don't have integer amount of cycles of your sinusoid. Please search for spectral leakage. Easy way to check that is to use the equation for frequency spacing in frequency domain, that is:

$$\Delta f=\dfrac{f_s}{N} $$

Where $f_s $ is your sampling frequency. Then your frequency vector is:

$$f=[0, \; \Delta f, \; 2\Delta f, \; ... \;,(N-2)\Delta f, \;(N-1)\Delta f] $$

If frequency you are analysing matches exactly one in your frequency vector, then amplitude will be correctly estimated, otherwise - leakage. In order to minimize this effect, you might want to apply window to your signal. It will reduce amount of energy leaking to other bins and in slight way will improve estimation of original amplitude, although you won't be really able to retrieve it in 100%.

Here is small modification of your code to help you visualize that. The red curve is ideal spectrum and blue one is for given number of frequency bins (strictly connected to your signal length).

import numpy as np
import matplotlib.pyplot as plt
import time

def get_fft(y, fs):
    """ Get the FFT of a given signal and corresponding frequency bins.

    Parameters:
        y  - signal
        fs - sampling frequency
    Returns:
        (mag, freq) - tuple of spectrum magitude and corresponding frequencies
    """
    n  = len(y)      # Get the signal length
    dt = 1/float(fs) # Get time resolution

    fft_output = np.fft.rfft(y)     # Perform real fft
    rfreqs = np.fft.rfftfreq(n, dt) # Calculatel frequency bins
    fft_mag = np.abs(fft_output)    # Take only magnitude of spectrum

    # Normalize the amplitude by number of bins and multiply by 2
    # because we removed second half of spectrum above the Nyqist frequency 
    # and energy must be preserved
    fft_mag = fft_mag * 2 / n           

    return np.array(fft_mag), np.array(rfreqs)

def build_signal(A, f, N, fs):
    """ Generate cosinusoidal signal of a given frequency and amplitude.

    Parameters:
        A   - amplitude of a signal
        f   - fundamental frequency
        N   - number of signal samples
        fs  - sampling frequency
    Returns:
        t   - time vector
        y   - generated signal    
    """

    dt = 1/float(fs)            # Get the time domain spacing
    T  = N*dt               # Upper time limit
    t  = np.arange(0, T, dt)    # Construct time vector

    y = A * np.cos(2*np.pi*f*t) # Generate the signal

    return t, y

if __name__ == "__main__":

    # Define the parameters
    fs = 100   # Sampling frequency (must satisfy the Nyquist theorem)
    f  = 10    # Signal fundamental frequency
    A  = 0.5   # Signal amplitude

    # Prepare figure for plotting
    plt.figure()
    plt.ion()
    plt.show(block=False)

    f_off = 1                     # Frequency offset in Hz for zoom

    # Precomute the spectrum without any leakage
    (t, y) = build_signal(A, f, 1000, fs)
    (MAG, FREQ) = get_fft(y, fs)

    for N in np.arange(1000, 10010, 1):

        # Clear and scale
        plt.clf()

        # Generate the signal  
        (t, y) = build_signal(A, f, N, fs)

        # Get the spectrum and frequencies
        (mag, freq) = get_fft(y, fs)

        # Plot the reference spectrum
        plt.plot(FREQ, MAG, '-or')
        plt.hold('on')
        # Plot the spectrum for given number of bins
        plt.plot(freq, mag, '-o')
        plt.title('Frequency spectrum for %d bins' % N)
        plt.xlabel('Frequency')
        plt.ylabel('Amplitude')
        plt.grid('on')
        plt.hold('off')
        plt.legend(['Ideal spectrum', 'Spectrum for %d bins' % N])
        plt.xlim((f-f_off, f+f_off))
        plt.ylim((0, A))      

        plt.draw()
        plt.pause(0.5)

Two intermediate steps. Signal is 10 Hz sinusoid with an amplitude of 0.5, sampling frequency is 100 Hz.

Exact match:

enter image description here

Frequency component you are analysis falls between two frequency bins:

enter image description here

$\endgroup$
  • $\begingroup$ I've combined two sine waves one having amplitude and frequency of 1000 and 390Hz and the second with 2000 and 440Hz respectively. I've used a sampling rate of 1000Hz but I can't retrieve the correct amplitude even after increasing bins. The most I could get is an amplitude of 1400 at 440Hz and 700 at 390Hz approximately. Is there a certain threshold after which leakage can't be avoided or am I doing something wrong? $\endgroup$ – Furqan Hashim Nov 17 '18 at 18:25
0
$\begingroup$

The maximum magnitude has to be interpolated for any sinusoid that is not exactly periodic by an exact integer number of periods in the FFT length. An interpolation is necessary because those frequencies are not at the center of any FFT result bins (integer multiples of Fs/N in frequency), but slightly between them, thus spreading out the peak energy among multiple result bins.

Just plotting the maximum magnitude of the FFT result vector won't give you that interpolation between bins. Thus you only see the magnitude of the nearest bin, which could be half a bin away from the real frequency peak, and thus with a smaller value than the sinusoid magnitude.

$\endgroup$
  • $\begingroup$ Could you provide me some hints about how I could make this interpolation ? $\endgroup$ – HadiM May 22 '14 at 20:13
  • $\begingroup$ Sinc reconstruction is the best interpolation kernel. See: ccrma.stanford.edu/~jos/resample/… $\endgroup$ – hotpaw2 Jul 11 '16 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.