1
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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()

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Note that you may need to interpolate the dominant frequency if it is between FFT result bins. If the dominant frequency is not exactly integer periodic in the FFT length, interpolating its phase (with reference to the Center of the FFT window) might more easily be done after an FFT shift. $\endgroup$ – hotpaw2 Oct 23 '14 at 1:23

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.