傅里叶变换（python）（2）例题

目的：去除噪音。时间和信号的关系是：gaussNoisy(t) （g(t)）

noiseAmp：控制噪音的大小。noiseAmp为0时图像会变成光滑的信号。

import numpy as np

def gaussNoisy(t, noiseAmp):
    noise = noiseAmp*(np.random.randn(len(t)))
    return np.exp(-0.5*t*t) * (1.0+noise)

N = 256
t = np.linspace(-4.0, 4.0, N)
g = gaussNoisy(t, 0.1)

（閾値は試行錯誤により決定しなさい：范围根据错误决定，我试了好多次xlim和ylim找噪音）

以下为解决方法的程序：

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as scifft
plt.close('all')
def gaussNoisy(t, noiseAmp):
    noise = noiseAmp*(np.random.randn(len(t)))
    return np.exp(-0.5*t*t) * (1.0+noise)

N = 256
t = np.linspace(-4.0, 4.0, N)
g = gaussNoisy(t, 0.1)


fig1, ax1 = plt.subplots(2, 1, figsize=(9, 7))
ax1[0].plot(t,g, color='red') 
ax1[0].set_xlim(-4, 4) 
ax1[0].set_ylim(0,1.3) 
ax1[0].set_xlabel('t') 
ax1[0].set_ylabel('g(t)')


num_data =len(t)
freq = scifft.fftfreq(num_data, d=0.03125)  
freq = scifft.fftshift(freq)

f_hat = scifft.fft(g,n=num_data)
f_hat = scifft.fftshift(f_hat) 
f_spect = f_hat * np.conj(f_hat)/num_data
fig1, ax2 = plt.subplots(4, 1, figsize=(9, 7))
ax2[0].plot(freq, np.real(f_hat), color='blue')
ax2[0].set_xlim(-128,128)
ax2[0].set_ylim(-1.2, 1.2) 
ax2[0].set_xlabel('freq (Hz)') 
ax2[0].set_ylabel('Re[f_hat]') 

# f_hat(f)作图（虚部imag）
ax2[1].plot(freq, np.imag(f_hat), color='blue')
ax2[1].set_xlim(-128,128) 
ax2[1].set_ylim(-1.2, 1.2) 
ax2[1].set_xlabel('freq (Hz)') 
ax2[1].set_ylabel('Im[f_hat]')
 
# f_hat的spect作图
ax2[2].plot(freq,f_spect, color='blue',label="before filtering") #傅里叶变换：可以看出有噪音
ax2[2].set_xlim(-128,128) 
ax2[2].set_ylim(-0.01, 0.01) 
ax2[2].set_xlabel('freq (Hz)') 
ax2[2].set_ylabel('f_spect') 
ax2[2].legend()

indices = f_spect>0.05 
#print("indices=", indices)
 

f_spect_filtered = f_spect * indices
f_hat_filtered = f_hat * indices
 
ax2[3].plot(freq,f_spect_filtered, color='red',label="after filtering")
ax2[3].set_xlim(-128,128)
ax2[3].set_ylim(-0.05, 0.05) 
ax2[3].set_xlabel('freq (Hz)') 
ax2[3].set_ylabel('f_spect(freq)') 
ax2[3].legend()
 
#逆傅里叶变换准备：要用fftshift
f_hat_filtered = scifft.fftshift(f_hat_filtered)
f_filterd = scifft.ifft(f_hat_filtered,n=num_data)
 
ax1[1].set_xlim(-4, 4) 
ax1[1].set_ylim(0,1.3) 
ax1[1].plot(t, np.real(f_filterd), color='blue') 
ax1[1].set_xlabel('time (sec)') 
ax1[1].set_ylabel('f_filtered(t)') 
 
plt.show()
plt.tight_layout()