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

发布于 2023-09-13  152 次阅读


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中,滤波后的频率-幅值图如下:

使用逆傅里叶变换可以得到处理后的信号:

注意上图中的第三幅图,将高通和低通滤波的结果进行逆变换、叠加之后获得了原始的有噪声的数据。

Everything not saved will be lost.
最后更新于 2023-09-13