Matlab FIR及IIR 数字滤波器 (一)

前言:

本文为本人在学习华科"数字信号分析理论与实践"课程所做的记录, 由于本人能力有限,在实现方法等各方面存在问题,欢迎交流指出.

本文将介绍实现最简单的FIR数字滤波器(fir1函数) ,以及附带的滤波函数(filter)

FIR 滤波器

原理分析

设计一个滤波器相当于设计一个系统,通过它的传递函数屏蔽掉不需要的输入信号,保留目标信息.

​ FIR(Finite Impulse Response)滤波器,即"有限长单位冲激响应滤波器",其脉冲响应是有限的持续时间，并在有限的时间内稳定为零。

其傅立叶系数h(n)实际上就是数字滤波器的冲击响应.

$h(k)=a(k),k=0,1,2,...,N$

其单边Z变换为

$H(z)=\sum^N_{k=0}{a_kZ^{-k}}$

• 低通滤波器

理想的滤波器频率响应方程为1

$H_d(\omega)=\begin{cases} 1& |\omega|\leq \omega_c \\ 0& \omega<|\omega|<\pi \end{cases}$

$h_d(k)=\frac{\omega_c}{pi}\frac{sin(\omega_ck)}{\omega_ck};k=0;\pm1,\pm2,..,\infty$

对系数进行归一化, 即

$\begin{aligned} a_0 &= v_1 = {f_c}/{\frac{1}{2}f_s}\\ a_n &= \frac{sin(n \pi v_1)}{n \pi};\ n=0;\pm1,\pm2,..,\infty \end{aligned}$

• 高通滤波器

理想的滤波器频率响应方程, 及系数计算

$\begin{aligned} &H_d(\omega)=\begin{cases} 1& |\omega|\geq \omega_c \\ 0& other \end{cases}\\ a_0 &= 1-v_1 \\ a_n &= -\frac{sin(n \pi v_1)} {n \pi};\ n=0;\pm1,\pm2,..,\infty \end{aligned}$

• 带通系数

理想的滤波器频率响应方程,及系数计算

$\begin{aligned} &H_d(\omega)=\begin{cases} 1& \omega\leq|\omega|\leq\pi \\ 0&other \end{cases}\\ a_0 &= v_2-v_1 \\ a_n &= \frac{sin(n \pi v_2)-sin(n \pi v_1)}{n \pi};\ n=0;\pm1,\pm2,..,\infty \end{aligned}$

由于h(k)可能是无限长序列，是物理不可实现的,为此要寻找一个h(n)，在相应的误差准则下最近逼近h(k)。

• 带阻系数

理想的滤波器频率响应方程,及系数计算

$\begin{aligned} &H_d(\omega)=\begin{cases} 0& \omega\leq|\omega|\leq\pi \\ 1&other \end{cases}\\ a_0 &= 1-(v_2-v_1) \\ a_n &= -\frac{sin(n \pi v_2)-sin(n \pi v_1)}{n \pi};\ n=0;\pm1,\pm2,..,\infty \end{aligned}$

• 系数处理(窗函数法)
  • $h(k)$ 可能是无限长序列，是物理不可实现的,为此要寻找一个 $h(n)$，在相应的误差准则下最近逼近$h(k)$。
  • 直接截取前n各系数, 故需要通过窗函数方法对系数进行平滑处理, 减小由于截断带来的跳动.

$h(n)=w(n)·a(n)$

本实现中默认使用汉宁(Hamming)窗函数对系数进行平滑处理.

程序实现

首先进行进行输入参数的判断

• 必须要输入的参数有采样点数目Np ,相对截止频率 v .其中v可以为二维数组
• 可选参数有 加窗类型 (通过指针p1保存) , 滤波器类型(通过指针p2保存)

function b = firl_fun(Np,v,varargin)

% % 使用方法 
% % firl_fun(Np,v) N为取样点数目，低通滤波器，截止频率v*fs
% % firl_fun(Np,v,@wintype) 加窗类型wintype,如 @hamming
% % firl_fun(Np,v,type) 滤波方式type可选'high'高通,'low'低通
% % firl_fun(Np,v,@wintype,type)
% % firl_fun(Np,[v1,v2]) 
% % firl_fun(Np,[v1,v2],type) 可选'bandpass'带通,'bandstop'带阻
% % firl_fun(Np,[v1,v2],@wintype,type)

N = floor(Np/2);
narginchk(2,4);%检查输入参数个数
p1 = @hamming; %p1指向窗函数(默认汉宁)
p2 = 'low'; %p2指向滤波类型（默认低通）

%判断滤波类型
switch nargin
    %只有3个参数的判断，
    case 3
        if isa(varargin{1},'function_handle')
            p1 = varargin{1};
        else
            p2 = varargin{1};          
        end
    %4参数读入    
    case 4
        p1 = varargin{1};
        p2 = varargin{2}; 
end

然后根据所选类型进行原理中的系数计算, 完成滤波器设计.

k=1; 
%判断滤波器类型
switch p2 
    case 'low'  %低通
        for i=-N:N 
            f=i*pi; 
            h(k)=sin(v(1)*f)/f; %计算响应函数
            k=k+1; 
        end
        h(N+1)=v(1);
    case 'high' %高通
        for i=-N:N 
            f=i*pi; 
            h(k)=-sin(v(1)*f)/f; 
            k=k+1; 
        end
        h(N+1)=1-v(1);
    case 'bandpass' %带通
        for i=-N:N 
            f=i*pi; 
            h(k)=(sin(v(2)*f)-sin(v(1)*f))/f; 
            k=k+1; 
        end  
        h(N+1)=v(2)-v(1);
    case 'bandstop' %带通阻止
        for i=-N:N 
            f=i*pi; 
            h(k)=-(sin(v(2)*f)-sin(v(1)*f))/f; 
            k=k+1; 
        end  
        h(N+1)=1-(v(2)-v(1));
end
b = h.*p1(2*N+1)';% -N:N的长度为2N+1
end

实际效果

低通、高通

带通、带阻

#### 测试代码

clear;
%设计的样本 %原信号
Fs = 2048;
dt=1.0/Fs;
T =1; N=T/dt;
t=[0:N-1]/N;

x1 =sin(2*pi*50*t)+sin(2*pi*300*t) +sin(2*pi*500*t);
subplot(3,2,1); 
plot(t,x1); 
axis([0, 0.1, -2,2]);
title("原信号时域"); 

P=fft(x1,N);
Pyy =2*sqrt(P.* conj(P))/N;
f=linspace(0,Fs/2,N/2);
subplot(3,2,2);
plot(f,Pyy(1:N/2));
title("原信号频域");

%选择滤波器
type = 4;
%--------------------------系统自带fir1 和filter
switch type 
    case 1 
        b = fir1(48,0.1);
    case 2 
        b = fir1(48,0.1,'high');
    case 3
        b = fir1(48,[0.1,0.3]);
    case 4
        b = fir1(48,[0.1,0.3],'stop');
end
x2= filter(b,1,x1);

%------------------------自己写的fir1_fun 和 filter_fun
switch type 
    case 1 
        b2 = firl_fun(48,0.1,@hamming);
    case 2 
        b2 = firl_fun(48,0.1,@hamming,'high');
    case 3
        b2 = firl_fun(48,[0.1,0.3],@hamming,'bandpass');
    case 4
        b2 = firl_fun(48,[0.1,0.3],@hamming,'bandstop');
end
x3= filter_fun(b2,1,x1);

%--------------------------（自带滤波器频域图）
subplot(3,2,3);
[h,w] = freqz(b,1,512);
plot(w/pi,abs(h),'linewidth',1)
title("Matlab的滤波器频域图"); 
%--------------------------（自带滤波器处理后）
subplot(3,2,5);
plot(t,x2);
axis([0, 0.1, -2,2]);
title("Matlab的滤波器处理后"); 

%--------------------------（编写的滤波器频域图）
subplot(3,2,4);
[h,w] = freqz(b2,1,512);
plot(w/pi,abs(h),'linewidth',1)
title("自写的滤波器频域图"); 
%--------------------------（编写的滤波器处理后）
subplot(3,2,6);
plot(t,x3);
axis([0, 0.1, -2,2]);
title("自写的滤波器处理后");

---

滤波器的实现

原理分析

按照Matlab的原理解释

那实现就很简单了

function y= filter_fun(b,a,f)
%这个滤波实际上就是把现有数据当作输入，滤波公式当作系统传递函数
%FIR时，利用卷积求出系统响应就好了，所及这个就相当于一个卷积算法。
%FII时，利用就根据差分方程计算出来。
n = length(b);
m = length(f);
na = length(a)-1;%分母最高阶次数
S = zeros(1,m);%卷积结果存放

%t=1:n-1 时为开头过度段，缺位补零相乘，求和时相当于没有该项
%t=n:m-1 时，按照卷积公式滑动求和
%t=m:m+n-1 时，时为尾部过度段，与 t=2:n 类似。
for t=1:n-1
    S(t) = sum(f(1:t).*fliplr(b(1:t))) ;
    if na >2
        if t > na
            S(t) = S(t) - sum( fliplr(a(2:na+1)).*S(t-1:t-na));
        else
            S(t) = S(t) - sum( fliplr(a(2:t)).*S(1:t-1));
        end
    end
end
for t = n:m-1
    S(t) = sum(f(t-n+1:t).*fliplr(b(1:n)));
    if na >2
        if t > na
            S(t) = S(t) - sum( fliplr(a(2:na+1)).*S(t-na:t-1));
        else
            S(t) = S(t) - sum( fliplr(a(2:t)).*S(1:t-1));
        end
    end    
    
end
y = S;