窗函数法设计FIR滤波器与IIR滤波器Python编写(从MATLAB移植）

窗函数法设计FIR滤波器与IIR滤波器Python编写(从MATLAB移植）

1. 题目要求

2. MATLAB代码
来自Peter_831

clc;
clear;
close all;

Fs=4000; % 采样频率4000Hz  
t=0:1/Fs:1;  
s50=sin(2*pi*50*t); % 产生50Hz正弦波  
s1000=sin(2*pi*1000*t); % 产生1000Hz正弦波  
s=s50+s1000; % 混合信号：信号叠加
  
figure; % 画图  
subplot(2,1,1);plot(s);grid on;  
axis([0 1000 -2 2]);
title('原始信号50Hz和1000Hz正弦波混合');  


% -----------------FFT分析信号频谱----------------------  
len = 1024;    
y=fft(s,len);  % 对原始输入信号做len点FFT变换  
f=Fs*(0:len/2 - 1)/len;  
  
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid on;  
title('原始信号频谱')  
xlabel('Hz');ylabel('幅值'); 


% ------------------------IIR模拟低通滤波器设计---------------------------  
Fp=200; % 通带截止频率200Hz  
Fc=1000; % 阻带截止频率1000Hz  
Rp=1; % 通带波纹最大衰减为1dB  
Rs=40; % 阻带衰减为40dB    20lg(1/100)=40dB   
%-----------------------计算最小滤波器阶数-----------------------------  
na=sqrt(10^(0.1*Rp)-1);  
ea=sqrt(10^(0.1*Rs)-1);  
N=ceil(log10(ea/na)/log10(Fc/Fp));   %巴特沃兹阶数  
%----------------------巴特沃兹低通滤波器------------------------------  
Wn=Fp*2/Fs;  
[Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器  
[BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线

figure;
A1 = BW*Fs/(2*pi);
plot(A1,abs(BH));
grid on;xlabel('频率／Hz'); ylabel('频率响应幅度'); 
title('巴特沃兹低通滤波器');


Bf=filter(Bb,Ba,s); % 进行低通滤波  
By=fft(Bf,len);  % 对滤波输出信号做len点FFT变换  

figure; % 画图  
subplot(2,1,1);plot(t*Fs,s50,'blue',t*Fs,Bf,'red');grid on;  
axis([0 1000 -1 1]);title('巴特沃兹低通滤波输出');
legend('50Hz正弦波信号','巴特沃斯滤波器滤波后');  
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid on;  
title('巴特沃斯低通滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  


%-----------------------确定相应的FIR数字滤波器指标-------------------------
wp=2*pi*1000/Fs;  %通带截止频率
% ws=2*pi*50/Fs;
ws=pi/4;      %阻带截止频率
Bt=wp-ws;     %过渡带带宽
N0=ceil(6.2*pi/Bt);  %汉宁窗窗长
N=N0+mod(N0+1,2);    %由于是高通滤波器，所以窗长N必须为奇数
wc=(wp+ws)/(2*pi);   %计算理想高通滤波器通带截止频率（关于pi归一化）
%----------------------FIR数字高通滤波器(汉宁窗)-------------------------
hn=fir1(N-1,wc,'high',hanning(N));  
[BH,BW]=freqz(hn,1); % 绘制频率响应曲线  

figure;
A1 = BW*Fs/(2*pi);
plot(A1,abs(BH));
grid on;xlabel('频率／Hz'); ylabel('频率响应幅度'); title('采用汉宁窗设置高通FIR数字滤波器');

Bf=filter(hn,1,s); % 进行高通滤波  
By=fft(Bf,len);  % 对高通滤波输出信号做len点FFT变换

figure; % 画图  
subplot(3,1,1);plot(t*Fs,s1000,'blue');grid on;  
axis([0 400 -1 1]);
title('原1000Hz正弦信号');
subplot(3,1,2);plot(t*Fs,Bf,'red');grid on;
axis([0 400 -1 1]);
title('FIR数字高通滤波后');
subplot(3,1,3);plot(f,abs(By(1:len/2)));grid;  
title('高通FIR数字滤波器滤波后信号频谱');  
xlabel('Hz');ylabel('幅值'); 

3. Python移植版本

from scipy.signal import butter
from scipy.signal import freqz
from scipy.signal import lfilter
from scipy.signal import hanning
from scipy.signal import firls
from scipy.signal import firwin
from matplotlib.pyplot import stem
import numpy as np
from numpy.fft import *

import matplotlib.pyplot as plt
import numpy.matlib 

from math import pi

from math import log2
from math import pi
from math import ceil
from math import sqrt
from cmath import exp
from math import log10

Fs=4000 # 采样频率4000Hz 
 
t=np.arange(0,1,1/Fs) 

c50=np.cos(2*pi*50*t) # 产生50Hz余弦波  
c1000=np.cos(2*pi*1000*t) # 产生1000Hz余弦波  
s=c50+c1000  # 信号叠加

#解决中文显示问题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
  
fig1=plt.figure(1)     
fig1.tight_layout() 

plt.subplot(2,1,1)
plt.plot(s)
 
plt.axis([0,1000,-2,2])
plt.grid()  
plt.title('原始信号50Hz和1000Hz余弦波混合')  


# -----------------FFT分析信号频谱----------------------  
lens = 1024    
y=fft(s,lens)  # 对原始输入信号做len点FFT变换  
f=Fs*np.arange(0,lens)/lens 
  
plt.subplot(2,1,2)
plt.plot(f,abs(y))
plt.grid()
plt.title('原始信号频谱')  
plt.xlabel('Hz')
plt.ylabel('幅值') 




# ------------------------IIR模拟低通滤波器设计---------------------------  
Fp=200 # 通带截止频率200Hz  
Fc=1000 # 阻带截止频率1000Hz  
Rp=1 # 通带波纹最大衰减为1dB  
Rs=40 # 阻带衰减为40dB  
#-----------------------计算最小滤波器阶数-----------------------------  
na=sqrt(10**(0.1*Rp)-1)  
ea=sqrt(10**(0.1*Rs)-1)  
N=ceil(log10(ea/na)/log10(Fc/Fp))   #巴特沃兹阶数  
#----------------------巴特沃兹低通滤波器------------------------------  
Wn=Fp*2/Fs  
[Bb,Ba]=butter(N,Wn,'low') 
[BW,BH]=freqz(Bb,Ba) 

fig2=plt.figure(2) 
A1 = BW*Fs/(pi)
plt.plot(A1,abs(BH))
plt.grid()
plt.xlabel('频率')
plt.ylabel('频响') 
plt.title('巴特沃兹低通滤波器')


Bf=lfilter(Bb,Ba,s) # 进行滤波  
By=fft(Bf,lens)  # 对滤波输出信号做FFT变换  

fig3=plt.figure(3)  # 画图  
plt.subplot(2,1,1)
plt.plot(t*Fs,c50,'blue',t*Fs,Bf,'red')
plt.grid 
plt.axis([0,1000,-1,1])
plt.title('巴特沃兹低通滤波输出')
plt.legend('50Hz余弦弦波信号','滤波后')  
plt.subplot(2,1,2)
plt.plot(f,abs(By))
plt.grid
plt.title('低通滤波后信号频谱')  
plt.xlabel('Hz')
plt.ylabel('幅值')  



#-----------------------确定相应的FIR数字滤波器指标-------------------------
wp=2*pi*1000/Fs  #通带截止频率
# ws=2*pi*50/Fs
ws=pi/4      #阻带截止频率
Bt=wp-ws     #过渡带带宽
N0=ceil(6.2*pi/Bt)  #汉宁窗窗长
N=N0+(N0+1)%2    #由于是高通滤波器，所以窗长N必须为奇数
wc=(wp+ws)/(2*pi)   #计算理想高通滤波器通带截止频率（关于pi归一化）
#----------------------FIR数字高通滤波器(汉宁窗)-------------------------
#hn=firls(N-1,wc,'high',hanning(N)) 
hn=firwin(N,wc,window='hamming',pass_zero='highpass') 

[BW,BH]=freqz(hn,1) # 绘制频率响应曲线  

fig4=plt.figure(4)  # 画图  
A1 = BW*Fs/(2*pi)
plt.plot(A1,abs(BH))
plt.grid()
plt.xlabel('频率') 
plt.ylabel('频响') 
plt.title('汉宁窗设置高通FIR数字滤波器')

Bf=lfilter(hn,1,s) # 进行高通滤波  
By=fft(Bf,lens)  # 对高通滤波输出信号做len点FFT变换

fig4=plt.figure(5)  # 画图 
fig4.tight_layout() 

plt.subplot(3,1,1)
plt.plot(t*Fs,c1000,'blue')
plt.grid 
plt.axis([0,400,-1,1])
plt.title('原余弦信号')
plt.subplot(3,1,2)
plt.plot(t*Fs,Bf,'red')
plt.grid
plt.axis([0,400,-1,1])
plt.title('高通滤波后')
plt.subplot(3,1,3)
plt.plot(f,abs(By))
plt.grid  
plt.title('高通FIR滤波频谱')  
plt.xlabel('频率')
plt.ylabel('幅值') 

4、运行结果




5、结果分析
1、Python的scipy库中有许多数字信号处理库，调用方式和MATLAB很像，不过有一些区别。要注意输出的参数相对应的次序。
比如：

%MATLAB
[BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线

#Python
[BW,BH]=freqz(Bb,Ba) 

2、解决matplotlib无法显示中文问题
加入以下语句

#解决中文显示问题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

3、Spyder的图片嵌入到编辑器中，分辨率低
解决办法：

记住设置完后记得重启。
重启、重启、重启，重要事情说三遍。