1 # -*- coding: utf-8 -*- 2 """ 3 Spyder Editor 4 5 This is a temporary script file. 6 """ 7 8 import matplotlib.pyplot as plt 9 import scipy as sci 10 import numpy as np 11 12 13 14 plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標簽 15 plt.rcParams['axes.unicode_minus']=False #用來正常顯示負號 16 17 18 fs = 200; #采樣頻率 19 f1 = 10; #信號分量1頻率 20 f2 = 20; #信號分量2頻率 21 f3 = 30; #信號分量3頻率 22 f4 = 40; #信號分量1頻率 23 n = np.arange(100)/fs; #采樣點 24 25 # 生成信號 26 yn = [20*np.sin(2*np.pi*f1*i) + 10*np.sin(2*np.pi*f2*i) + 5*np.sin(2*np.pi*f3*i) + 20*np.sin(2*np.pi*f4*i) for i in n]; 27 28 29 30 # fft分析原始信號 31 fft_result_abs1 = np.abs(sci.fft(yn))*2/len(n); 32 f_idx = np.arange(len(n))*fs/len(n); 33 34 # 濾波器設計和濾波 35 (b,a) = sci.signal.butter(N=5,Wn=(15),btype='lowpass',output='ba',fs=200); #生成低通濾波器參數 36 #print(b) 37 #print(a) 38 yn_filtered = sci.signal.filtfilt(b,a,yn); #濾波器 39 # fft分析濾波后的信號 40 fft_result_abs2 = np.abs(sci.fft(yn_filtered))*2/len(n); 41 42 43 44 # 畫圖 45 plt.close('all'); 46 plt.figure(); 47 plt.subplot(2,2,1); 48 plt.plot(n,yn); 49 plt.title(U'原始時域信號'); 50 plt.subplot(2,2,2); 51 plt.plot(f_idx[range(int(len(n)/2))],fft_result_abs1[range(int(len(n)/2))]); 52 plt.title(U'原始信號頻域分析'); 53 54 plt.subplot(2,2,3); 55 plt.plot(n,yn_filtered); 56 plt.title(U'低通濾波后時域信號'); 57 58 59 plt.subplot(2,2,4); 60 plt.plot(f_idx[range(int(len(n)/2))],fft_result_abs2[range(int(len(n)/2))]); 61 plt.title(U'低通濾波信號頻域分析'); 62 63 64 plt.figure(); 65 (w,h) = sci.signal.freqs(b,a); 66 #plt.plot(w, 20 * np.log10(abs(h))) #繪制幅頻響應，頻率軸取對數，幅度軸轉換成dB。 67 68 plt.plot(w, abs(h))
圖片顯示：
Ref：
https://blog.csdn.net/qq7835144/article/details/88838537
https://segmentfault.com/a/1190000005144275