Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------
理论部分见博文： Fourier Transform:1D Transform
1. 生成数据集
使用正弦三角函数生成一组2个周期、每个周期100个点（共200个点）的信号数据集，并向生成的数据中引入正态分布的噪声：
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 = fig1.add_subplot(111)
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. 无噪声数据集的处理
这里通过对无噪声数据集的处理来熟悉一维傅里叶变换过程和Scipy中各行代码的含义：
# 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 = fig2.add_subplot(111)
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()
将信号由time domain转换至frequency domain中后，对经傅里叶变换得到的复数的实部和虚部进行绘图如下：
在Frequency domain中的幅值、频率和相位角可以使用以下方式得到：
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
使用Scipy.fftpack.ifft()对信号进行逆变换就可以重新回到time domain:
# 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 = fig3.add_subplot(111)
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()
上图中黄色曲线是使用幅值（amplitude，缩小了100倍）、频率和相位角信息带入正弦函数得到的曲线，可以看到，在time domain 和 Frequency domain中的信号并不相同。
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()
ax4 = fig4.add_subplot(111)
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()
在Frequency domain中，滤波后的频率-幅值图如下：
使用逆傅里叶变换可以得到处理后的信号：
注意上图中的第三幅图，将高通和低通滤波的结果进行逆变换、叠加之后获得了原始的有噪声的数据。
Comments NOTHING