[Python] Fourier analysis

因为在UCLA Physics 180N上遇到了这样的project， 是对用python做傅立叶分析很好的练习，在这里分享一下以备用。

Fourier Transforms: 
One can always look at a physical process in time domain  or frequency domain . The powerful Fourier transform says that any periodic function  can be decomposed into sum of sine waves:

where the complex number  gives the amplitude and phase of the corresponding wave with . Equation  is the forward Fourier transforms and equation  is the backward Fourier transform.

If a signal  is discrete, sampled with  points, then we use discrete Fourier transforms (DFT):

For example: Given some data called sunspot, which contains monthly mean total sunspot number over the span of 50 years (so, it's discrete data points), we can use numpy.fft to perform fast Fourier transform:

sunspot_dft = np.fft.fft(sunspot)

Power spectra (power spectrum density) can be used when we want to see the distribution of the wave amplitudes over a range of frequency, which can reveal information not obvious in the time domain:

The peaks in the spectra are the dominant frequency modes that are most essential to describing the signal. So, a white noise signal will give a flat power spectrum. 

Here we can also take power spectrum of the sunspot data:

sunspot_ps = (T/N)**2 * (np.abs(sunspot_dft))**2 # power spetrum

The correlation of two signals X(t) and Y(t) is:

which tells us how well the two signals match when one is shifted by time lag . The correlation is a function if , so it lies in the time domain. When , the correlation function is known as autocorrelation function.

With DFT in equation  and , we can rewrite the correlation function as:

The signals can be normalized using:

#Normalize sunspot data 
sunspot_mean = np.mean(sunspot)
sunspot_std = np.std(sunspot)
sunspot_norm =(sunspot - sunspot_mean)/sunspot_std
sunspot_auto = np.correlate(sunspot_norm, sunspot_norm, mode='full' )
plt.plot(times, sunspot_auto/N)
plt.xlabel('Time lag (months)')
plt.ylabel('Autocorrelation')
plt.grid()
plt.show()

Lower pass filtering to clean up signal: 

Actually, the raw sunspot and cosmic ray count data is quite noisy, we can apply apply a low-pass filter to smooth out the signals. This is done by setting a cutoff frequency  and setting all signal frequencies higher than the cutoff to 0. We first transform the signal to frequency domain via Fourier transform, filter out the modes with frequency larger than cutoff frequency, and the filtered signal back the time domain via inverse Fourier transform

1. '''Low pass noise filter'''
2. def lowpass_filter(dft_data, cutoff):
3. filtered =np.zeros(N, dtype = np.complex_)
4. for n in np.arange(cutoff):
5. filtered[n] = dft_data[n]
6. return np.real(np.fft.ifft(filtered)) # ifft to transform back to time-domain

Here is what it looks like after apply low-pass filter of different cut-off frequencies: