赞
踩
最近从MATLAB转到python进行科学计算,来比较下MATLAB和SciPy计算巴特沃斯通带滤波的效果。
clear; close all; clc; sampRat = 100.0; %采样频率 T = 10; t = 0:1/sampRat:(T*sampRat-1)/sampRat; x = 5*sin(2*pi*t*5)+10*sin(2*pi*t*25); ff = fft(x); ff = abs(ff); ff = ff*2/T/sampRat; f = 0.0:1/T:(sampRat-1/T); figure() plot(f, ff); title('滤波前'); [b,a] = butter(4,[0.01/(sampRat/2),15/(sampRat/2)]);%设计butterworth滤波器(带通) y =filter(b,a,x); fff = fft(y); fff = abs(fff); fff = fff*2/T/sampRat; figure() plot(f, fff) title('滤波后');
SciPy
# -*- coding: utf-8 -*- import scipy.signal as signal import matplotlib.pyplot as plt import numpy as np import math from scipy.signal import butter, lfilter def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): b, a = butter_bandpass(lowcut, highcut, fs, order=order) y = lfilter(b, a, data) return y def butter_bandpass(lowcut, highcut, fs, order=5): nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(order, [low, high], btype='bandpass') return b, a sampRat = 100 dt = 1/sampRat T = 10 t = np.linspace(0, T, T*sampRat, endpoint=False) y = 5*np.sin(2*math.pi*t*5)+10*np.sin(2*math.pi*t*25) filterY = butter_bandpass_filter(y, 0.01, 15, sampRat, 4) f = np.linspace(0, sampRat, T*sampRat, endpoint=False) ff = np.fft.fft(y) ff = np.abs(ff)*2/T/sampRat plt.figure() plt.plot(f, ff) plt.title('滤波前的频谱') plt.show() ff = np.fft.fft(filterY) ff = np.abs(ff)*2/T/sampRat plt.figure() plt.plot(f, ff) plt.title('滤波后的频谱') plt.show()
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。