Find a sinusoidal fit from data with fft

Question

I have a vector with raw data (120 samples). From the frequency spectrum I can tell there is one dominant frequency component. Now I am supposed to use the fft data (magnitude and phase) to find a good sinusoidal fit and offset to fit the original data.

I am completely stuck on this one. Can anyone help?

Answer 1

Depending on kind of signal you are dealing with I can think of the easiest method to do so.

1. Take the DFT of your signal (remember about proper scaling of values).
2. Find the dominant frequency $f_0$ and its amplitude $A$ at magnitude spectrum.
3. From phase spectrum find the corresponding phase shift $\theta$ for this given frequency - usually it is given in radians.
4. Reconstruct your signal as: $x(t)=A\sin(2\pi f_0t+\theta) $

Some quick and dirty code is included below. No warranty of working whatsoever.

#!/usr/bin/env python

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

if __name__ == "__main__":
    # Signal and analysis parameters
    fs = 4000.0    # Sampling frequency
    T  = 2         # Duration of a signal
    f0 = 100       # Fundamental frequency
    A  = 1         # Sine amplitude
    #th = np.pi/4   # Phase shift in radians (45 degrees)
    th = 0

    # Time vector
    t = np.arange(0, T*fs)/fs 

    # Generate the signal with higher harmonics
    x = A*np.sin(2*np.pi*f0*t + th) + \
        A*np.sin(2*np.pi*2*f0*t + 2*th)/8 + \
        A*np.sin(2*np.pi*3*f0*t + 3*th)/16

    # Add some noise 
    x += 0.5*np.random.randn(x.size)    

    # Number of all samples
    N = x.size    

    # Apply the window
    win = sig.gaussian(N, N)
    x = x*win    

    # Perform the DFT
    X = np.fft.fft(x)
    # Get the frequency vector
    freq = np.linspace(0, fs-1/fs, N)

    # Get the phase and magnitude
    Xmag = np.abs(X)/np.sum(win)*2
    Xph  = np.angle(X)

    # Magnitude in logarithmic scale
    Xmagl = 10*np.log10(Xmag/Xmag.max())

    # Find the amplitude of the dominant frequency    
    ind_max = np.argmax(Xmag[0:np.ceil(N/2)])

    # Extract the amplitude ...
    A_ = Xmag[ind_max]
    # ... frequency ..
    f0_ = freq[ind_max]
    # ... phase
    th_ = Xph[ind_max] + np.pi/2

    # Reconstruct the signal
    x_ = A_ * np.sin(2*np.pi*f0_*t + th_)

    # Do some plotting
    plt.subplot(311)
    plt.plot(t, x)
    plt.grid(True)  
    plt.hold(True)
    plt.plot(t, x_, 'r')
    plt.title('Time domain signal')

    plt.subplot(312)
    plt.plot(freq[0:N/2], Xmagl[0:N/2])
    plt.grid(True)
    plt.title('Magnitude')

    plt.subplot(313)
    plt.plot(freq[0:N/2], Xph[0:N/2])
    plt.grid(True)
    plt.title('Phase')
    plt.show()