MATLAB和scipy计算巴特沃斯通带滤波的比较

引言

最近从MATLAB转到python进行科学计算，来比较下MATLAB和SciPy计算巴特沃斯通带滤波的效果。

MATLAB

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()