低通滤波器实现过程解析

首先来看一个低通滤波的例子：给定一组波形如下：

      如果想要获取低通0-0.5Hz的波形：需要调用低通滤波器可得到：

很简单的一个例子：程序见：点击打开链接

低通滤波器实现过程：

主函数：

y=lowp(x,f1,f3,rp,rs,Fs)

X:输入信号；Fs: 采样频率 f1:通带截止频率，rp:通带衰减 f3:阻带截止频率，rs:阻带衰减

Step1：模拟频率/数字频率转换

wp=2*pi*f1/Fs; ws=2*pi*f3/Fs;

 

Step2：求Chebyshev Type I filter order

[n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs); N:滤波器阶数 Wn：滤波器截止频率

wp/pi:归一化通带截止频率 ws/pi：归一化阻带截止频率 rp：通带衰减 rs:阻带衰减

 

Step3：求滤波器分子、分母系数

[bz1,az1]=cheby1(n,rp,wp/pi); bz1：分子系数矩阵； az1：分母系数矩阵

N：阶数 Rp：通带衰减 wp/pi：归一化通带衰减频率

 

 

Step3.1：get analog, pre-warped frequencies

 

if ~analog,     fs = 2;     u = 2*fs*tan(pi*Wn/fs); else     u = Wn; end

 

Step3.2:convert to low-pass prototype estimate

if btype == 1   % lowpass     Wn = u; elseif btype == 2   % bandpass     Bw = u(2) - u(1);     Wn = sqrt(u(1)*u(2));   % center frequency elseif btype == 3   % highpass     Wn = u; elseif btype == 4   % bandstop     Bw = u(2) - u(1);     Wn = sqrt(u(1)*u(2));   % center frequency end

 

 

Step3.3:Get N-th order Chebyshev type-I lowpass analog prototype

 

[z,p,k] = cheb1ap(n, r); Z:滤波器零点 P:滤波器极点 K:系数	N:滤波器阶数 r:通带衰减	包含：asinh，exp，flipud，complex，prod
[a,b,c,d] = zp2ss(z,p,k); 状态方程     x = Ax + Bu              y = Cx + Du   A,b,c,d:状态方程的系数	z,p,k：同上	包含：parse_input， zp2tf，tf2ss，isfinite cplxpair，poly. 待处理：zp2tf，tf2ss 异常处理：Try...catch	零极点转化为状态方程

 

Step3.4:Transform to lowpass, bandpass, highpass, or bandstop of desired Wn

[a,b,c,d] = lp2lp(a,b,c,d,Wn);

a,b,c,d：同上 Wn:截止频率

Step3.5 ：Use Bilinear transformation to find discrete equivalent:

[a,b,c,d] = bilinear(a,b,c,d,fs);	双线性变换
den = poly(a);	分子
num = cheb1num(btype,n,Wn,Bw,analog,den,r);	分母

 

Step4：差分方程求滤波后的序列

 

y=filter(bz1,az1,x);

bz1：分子系数矩阵； az1：分母系数矩阵 X：输入序列

 -------------------------------------------------------------------------------------

matlab实现是相当简单的，问题难得是如果要用C语言实现肯定是要费一番功夫的。

under-edit