# Find a sinusoidal fit from data with fft

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?

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

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


• 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. – hotpaw2 Oct 23 '14 at 1:23