python时域信号处理Demo:fft psd 自相关分析

东方耀

python时域信号处理Demo:fft psd 自相关分析

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. 
   
4. # python时域信号处理Demo:fft psd 自相关分析
   
5. N = 200     # 采样个数
   
6. f_s = 50    # 采样频率
   
7. t = np.linspace(0.0, N/f_s, N)
   
8. print("4秒内多少个采样点：", len(t))
   
9. signal1 = 2 * np.sin(2 * 1 * np.pi * t)   # f=1
   
10. signal2 = np.sin(2 * 5 * np.pi * t)       # f=5
    
11. # 时域信号合成
    
12. signal = signal1 + signal2
    
13. 
    
14. 
    
15. def get_values(y_values, N, f_s):
    
16.     f_values = np.linspace(0.0, N/f_s, N)
    
17.     return f_values, y_values
    
18. 
    
19. 
    
20. 
    
21. plt.figure(figsize=(16, 12))
    
22. 
    
23. plt.subplot(231)
    
24. f_values, signal_values = get_values(signal1, N, f_s)
    
25. plt.plot(f_values, signal_values, linestyle='-', color='blue')
    
26. plt.xlabel('Time [sec]', fontsize=16)
    
27. # Amplitude 振幅
    
28. plt.ylabel('Amplitude', fontsize=16)
    
29. plt.title("Time domain of the signal1(f=1)", fontsize=16)
    
30. # plt.show()
    
31. 
    
32. plt.subplot(232)
    
33. f_values, signal_values = get_values(signal2, N, f_s)
    
34. plt.plot(f_values, signal_values, linestyle='-', color='blue')
    
35. plt.xlabel('Time [sec]', fontsize=16)
    
36. plt.ylabel('Amplitude', fontsize=16)
    
37. plt.title("Time domain of the signal2(f=5)", fontsize=16)
    
38. # plt.show()
    
39. 
    
40. plt.subplot(233)
    
41. f_values, signal_values = get_values(signal, N, f_s)
    
42. plt.plot(f_values, signal_values, linestyle='-', color='blue')
    
43. plt.xlabel('Time [sec]', fontsize=16)
    
44. plt.ylabel('Amplitude', fontsize=16)
    
45. plt.title("Time domain of the signal", fontsize=16)
    
46. # plt.show()
    
47. 
    
48. 
    
49. from scipy.fftpack import fft
    
50. 
    
51. def get_fft_values(y_values, N, f_s):
    
52.     # 信号频率不会大于采样频率的一半
    
53.     f_values = np.linspace(0.0, f_s / 2.0, N // 2)
    
54.     fft_values_ = fft(y_values)
    
55.     print("原始信号数据经过fft后：", len(fft_values_))
    
56.     # 对称性 求模（振幅谱） 归一化(为什么乘以2.0)
    
57.     fft_values = 2.0 / N * np.abs(fft_values_[0:N // 2])
    
58.     print("单边振幅谱(归一化) 为什么乘以2.0???")
    
59.     print("绝对值就是求模，归一化后的fft值：", len(fft_values), fft_values[:5])
    
60.     return f_values, fft_values
    
61. 
    
62. 
    
63. plt.subplot(234)
    
64. f_values, fft_values = get_fft_values(signal, N, f_s)
    
65. plt.plot(f_values, fft_values, linestyle='-', color='blue')
    
66. plt.xlabel('Frequency [Hz]', fontsize=16)
    
67. plt.ylabel('Amplitude', fontsize=16)
    
68. plt.title("The FFT of our composite signal", fontsize=16)
    
69. # plt.show()
    
70. 
    
71. 
    
72. from scipy.signal import welch
    
73. 
    
74. 
    
75. def get_psd_values(y_values, N, f_s):
    
76.     # 用韦尔奇方法估计功率谱密度
    
77.     f_values, psd_values = welch(y_values, fs=f_s)
    
78.     return f_values, psd_values
    
79. 
    
80. 
    
81. plt.subplot(235)
    
82. # PSD英文全称Power Spectral Density，即功率谱密度，它和FFT一样，反映的是信号在频域上的信息。
    
83. # 其中PSD频谱图脉冲下方的面积表示信号在该频率上的能量分布。
    
84. f_values, psd_values = get_psd_values(signal, N, f_s)
    
85. plt.plot(f_values, psd_values, linestyle='-', color='blue')
    
86. plt.xlabel('Frequency [Hz]', fontsize=16)
    
87. plt.ylabel('PSD [V**2 / Hz]', fontsize=16)
    
88. plt.title("The PSD of our composite signal", fontsize=16)
    
89. 
    
90. 
    
91. def autocorr(x):
    
92.     # numpy互相关函数np.correlate
    
93.     # 《Python 数字信号处理应用》这本书上刚好看到过，题主想要用的应该是np.corrcoef这个方法，np.corrlated，是信号处理领域的自相关算法，
    
94.     # 【使用的非相关的标准化版定义，随着时滞变大，两个信号之间重叠点数减少，所以相关幅度降低了】，这是书上原话，
    
95.     # 我的理解是由于np.corrlated在有限范围内计算信号自相关，也就是错位相乘，随着时延增加重叠的点就会减少，计算的自相关也会变小，
    
96.     # 相比np.corrcoef，自相关函数的幅值明显有下降，需要除以N/2的长度这样可以获得与np.corrcoef几乎相同的结果
    
97.     result = np.correlate(x, x, mode='full')
    
98.     return result[len(result) // 2:]
    
99. 
    
100. 
     
101. def get_autocorr_values(y_values, N, f_s):
     
102.     autocorr_values = autocorr(y_values)
     
103.     # 细节信息
     
104.     x_values = np.array([1.0 * jj / f_s for jj in range(0, N)])
     
105.     return x_values, autocorr_values
     
106. 
     
107. 
     
108. plt.subplot(236)
     
109. # Autocorrelation是自相关的意思，它可以求出信号的自相关性，即信号经过一个时延后与自身的相似性
     
110. # 有趣的是Autocorrelation与PSD是一组FFT变换对，对Autocorrelation作FFT变换可得到PSD，对PSD作IFFT（快速傅里叶逆变换）可得到Autocorrelation
     
111. t_values, autocorr_values = get_autocorr_values(signal, N, f_s)
     
112. plt.plot(t_values, autocorr_values, linestyle='-', color='blue')
     
113. plt.xlabel('time delay [s]', fontsize=16)
     
114. plt.ylabel('Autocorrelation amplitude', fontsize=16)
     
115. plt.title("The autocorrelation of signal.", fontsize=16)
     
116. 
     
117. # plt.savefig('fft_psd_corr.png', dpi=200)
     
118. plt.show()
     

复制代码