# [Scipy] 使用Scipy探索一维傅里叶变换(Fourier Transform 1D)

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

## 1. 生成数据集

from scipy import fftpack
import numpy as np
import matplotlib.pyplot as plt

"""Make up some random data"""
time_step = 0.01
time_vec = np.arange(0, 2, time_step)

# A signal with a random normal distributed noise (2 cycles)
sig = np.sin(2 * np.pi * time_vec)
np.random.seed(1234)
sig_noise = sig + 0.25 * np.random.randn(time_vec.size)

# plot the signal before and after adding noise
fig1 = plt.figure()
ax1.plot(time_vec, sig, label='Original signal')
ax1.plot(time_vec, sig_noise, label='Noisy signal')
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Amplitude')
ax1.legend(loc='best')
ax1.plot()
plt.show()


## 2. 无噪声数据集的处理

# FFT of the signal without noise
sig_fft = fftpack.fft(sig)
# real part of the complex number
real_sig_fft = sig_fft.real
# imaginary part of the complex number
imag_sig_fft = sig_fft.imag

fig2 = plt.figure()
ax2.plot(real_sig_fft, imag_sig_fft, label='Original signal')
# ax2.plot(real_sig_noise_fft, imag_sig_noise_fft, label='Noisy signal')
ax2.set_xlabel('Real part')
ax2.set_ylabel('Imaginary part')
ax2.legend(loc='best')
plt.show()


amplitude = np.abs(sig_fft)  # amplitude spectrum (magnitude)
f = fftpack.fftfreq(sig.size, d=time_step)  # return the Discrete Fourier Transform sample frequencies
phase_angle = np.angle(sig_fft)  # return the phase angle of the signal


# Find the frequency, amplitude and phase angle of the signal with maximum amplitude
amp_peak_index = amplitude.argmax()  # index of the maximum amplitude
amp_peak = amplitude[amp_peak_index]  # maximum amplitude
freq_peak = f[amp_peak_index]  # frequency of the maximum amplitude
phase_angle_peak = phase_angle[amp_peak_index]  # phase angle of the maximum amplitude


# apply inverse FFT to the signal without noise
sig_ifft = fftpack.ifft(sig_fft)  # inverse FFT

# plot a signal with amplitude, phase and frequency
fig3 = plt.figure()
ax3.plot(time_vec, sig, label='Original signal')
ax3.plot(time_vec, amp_peak/100 * np.sin(2 * np.pi * freq_peak * time_vec + phase_angle_peak), label='Frequency domain (Amplitude/100)')
ax3.plot(time_vec, sig_ifft, label='Inverse FFT', linestyle=':', color='red')
ax3.set_xlabel('Time [s]')
ax3.set_ylabel('Amplitude')
ax3.legend(loc='best')
plt.show()


## 3. 有噪声数据集的处理

# denoise the noisy signal by removing high frequencies
sig_noise_fft = fftpack.fft(sig_noise)
amplitude_noise = np.abs(sig_noise_fft)
f_noise = fftpack.fftfreq(sig_noise.size, d=time_step)
phase_angle_noise = np.angle(sig_noise_fft)

# Find the frequency, amplitude and phase angle of the signal with maximum amplitude
amp_peak_index_noise = amplitude_noise.argmax()  # index of the maximum amplitude
amp_peak_noise = amplitude_noise[amp_peak_index_noise]  # maximum amplitude
freq_peak_noise = f_noise[amp_peak_index_noise]  # frequency of the maximum amplitude
phase_angle_peak_noise = phase_angle_noise[amp_peak_index_noise]  # phase angle of the maximum amplitude

# High pass filter
cutoff_idx = amplitude_noise < (amplitude_noise.max() / 6)
sig_noise_high_pass = sig_noise_fft.copy()
sig_noise_high_pass[cutoff_idx] = 0
f_noise_high_pass = f_noise[cutoff_idx]
amplitude_noise_high_pass = amplitude_noise[cutoff_idx]

# Low pass filter
cutoff_idx = amplitude_noise > (amplitude_noise.max() / 6)
sig_noise_low_pass = sig_noise_fft.copy()
sig_noise_low_pass[cutoff_idx] = 0
f_noise_low_pass = f_noise[cutoff_idx]
amplitude_noise_low_pass = amplitude_noise[cutoff_idx]

fig4 = plt.figure()
ax.set_title('FFT filter')
ax4.plot(f_noise, amplitude_noise, label='Noisy signal')
ax4.plot(f_noise_high_pass, amplitude_noise_high_pass, label='High pass filter')
ax4.plot(f_noise_low_pass, amplitude_noise_low_pass, label='Low pass filter')
ax4.set_xlabel('Frequency [Hz]')
ax4.set_ylabel('Amplitude')
ax4.legend(loc='best')
plt.show()


Everything not saved will be lost.