CODEWITHSUNDEEP
pythonnumpyfftcurve-fitting
andrerud
andrerud

Reputation: 166

Fourier series data fit with numpy: fft vs coding

Suppose I have some data, y, to which I would like to fit a Fourier series. On this post, a solution was posted by Mermoz using the complex format of the series and "calculating the coefficient with a riemann sum". On this other post, the series is obtained through the FFT and an example is written down.

I tried implementing both approaches (image and code below - notice everytime the code is run, different data will be generated due to the use of numpy.random.normal) but I wonder why I am getting different results - the Riemann approach seems "wrongly shifted" while the FFT approach seems "squeezed". I am also not sure about my definition of the period "tau" for the series. I appreciate the attention.

I am using Spyder with Python 3.7.1 on Windows 7

Example

import matplotlib.pyplot as plt
import numpy as np

# Assume x (independent variable) and y are the data.
# Arbitrary numerical values for question purposes:
start = 0
stop = 4
mean = 1
sigma = 2
N = 200
terms = 30 # number of terms for the Fourier series

x = np.linspace(start,stop,N,endpoint=True) 
y = np.random.normal(mean, sigma, len(x))

# Fourier series
tau = (max(x)-min(x)) # assume that signal length = 1 period (tau)

# From ref 1
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*x/tau)
    return c.sum()/c.size
def f(x, Nh):
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(1,Nh+1)])
    return f.sum()
y_Fourier_1 = np.array([f(t,terms).real for t in x])

# From ref 2
Y = np.fft.fft(y)
np.put(Y, range(terms+1, len(y)), 0.0) # zero-ing coefficients above "terms"
y_Fourier_2 = np.fft.ifft(Y)

# Visualization
f, ax = plt.subplots()
ax.plot(x,y, color='lightblue', label = 'artificial data')
ax.plot(x, y_Fourier_1, label = ("'Riemann' series fit (%d terms)" % terms))
ax.plot(x,y_Fourier_2, label = ("'FFT' series fit (%d terms)" % terms))
ax.grid(True, color='dimgray', linestyle='--', linewidth=0.5)
ax.set_axisbelow(True)
ax.set_ylabel('y')
ax.set_xlabel('x')
ax.legend()

Upvotes: 4

Views: 4707

Answers (1)

francis
francis

Reputation: 9817

Performing two small modifications is sufficient to make the sums nearly similar to the output of np.fft. The FFTW library indeed computes these sums.

1) The average of the signal, c[0] is to be accounted for:

f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(0,Nh+1)]) # here : 0, not 1

2) The output must be scaled. 

y_Fourier_1=y_Fourier_1*0.5

enter image description here The output seems "squeezed" because the high frequency components have been filtered. Indeed, the high frequency oscillations of the input have been cleared and the output looks like a moving average.

Here, tau is actually defined as stop-start: it corresponds to the length of the frame. It is the expected period of the signal.

If the frame does not correspond to a period of the signal, you can guess its period by convoluting the signal with itself and finding the first maximum. See Find period of a signal out of the FFT Nevertheless, it is unlikely to work properly with a dataset generated by numpy.random.normal : this is an Additive White Gaussian Noise. As it features a constant power spectral density, it can hardly be discribed as periodic!

Upvotes: 3

Related Questions