这篇文章本意就是做出这样的一个过程:
一个被50Hz噪声所干扰的10Hz信号,在时域很难被去躁并分离出来,但如果我们用DFT将其转换为频域的信号,就可以很容易的使用滤波器完成分离了。
那么上述过程中的四张图片,true signal,noise signal,input signal和频域内的input signal,我们如何通过python作出呢?下面给出代码及生成的信号图:
import numpy as np import matplotlib.pyplot as plt %matplotlib inline def generate_sinusoid(N, A, f0, fs, phi): ''' N(int) : number of samples A(float) : amplitude f0(float): frequency in Hz fs(float): sample rate phi(float): initial phase return x (numpy array): sinusoid signal which lenght is M ''' T = 1/fs n = np.arange(N) # [0,1,..., N-1] x = A * np.sin( 2*f0*np.pi*n*T + phi ) return x def generate_sinusoid_2(t, A, f0, fs, phi): ''' t (float) : time length of the generated sequence A (float) : amplitude f0 (float) : frequency fs (float) : sample rate phi(float) : initial phase returns x (numpy array): sinusoid signal sequence ''' T = 1.0/fs N = t / T return generate_sinusoid(N, A, f0, fs, phi) A = 2 f0 = 10 fs = 500 phi = 0 t = 1 x = generate_sinusoid_2(t, A, f0, fs, phi) n = np.arange(0, t, 1/fs) plt.plot(n, x) plt.ylim(-2.5, 2.5) plt.xlabel('time(s)') plt.ylabel('amplitude') plt.title('True signal in time domain') %config InlineBackend.figure_format = 'svg' plt.savefig('e_DFT1.eps', bbox_inches='tight') plt.show() A = 0.3 f0 = 50 fs = 500 phi = 0 t = 1 x = generate_sinusoid_2(t, A, f0, fs, phi) n = np.arange(0, t, 1/fs) plt.plot(n, x) plt.ylim(-2.5, 2.5) plt.xlabel('time(s)') plt.ylabel('amplitude') plt.title('Noise signal in time domain') %config InlineBackend.figure_format = 'svg' plt.savefig('e_DFT2.eps', bbox_inches='tight') plt.show() A = 2 f0 = 10 fs = 500 phi = 0 t = 1 y = generate_sinusoid_2(t, A, f0, fs, phi) n = np.arange(0, t, 1/fs) plt.plot(n, x+y) plt.ylim(-2.5, 2.5) plt.xlabel('time(s)') plt.ylabel('amplitude') plt.title('Input signal in time domain') %config InlineBackend.figure_format = 'svg' plt.savefig('e_DFT3.eps', bbox_inches='tight') plt.show() from scipy.fftpack import fft # fft is X = fft(x+y) mX = np.abs(X) # magnitude pX = np.angle(X) # phase #plt.subplot(2,1,1) plt.plot(mX/250) plt.xlabel('frequency(Hz)') plt.ylabel('amplitude') plt.title('Input signal in frequency domain') %config InlineBackend.figure_format = 'svg' plt.savefig('e_DFT4.eps', bbox_inches='tight') #plt.subplot(2,1,2) #plt.plot(pX) plt.show()参考文章: 信号生成及DFT的python实现: https://blog.csdn.net/weiwei9363/article/details/83653295