6
25
2011
0

IIR数字滤波器的设计与滤波

 

设计IIR滤波器,实现对存在加性干扰的时域离散信号进行滤波。

已知带加性干扰的信号用x(n)表示,x(n)=xs(n)+η(n),式中xs(n)是有用的信号,是一个0~0.2πrad的带限信号。η(n)是一个干扰信号,其频谱分布在0.3πrad以上。要求设计一个巴特沃斯IIR数字滤波器对信号x(n)进行滤波,将干扰η(n)滤除。要求在xs(n) 所在的通带内滤波器幅度平坦,在0.2πrad处幅度衰减不大于1dB,在噪声所在的0.3πrad 以上的频带内滤波器幅度衰减大于等于40 dB。

步骤:1.根据题目要求确定要设计的数字滤波器的技术指标(低通滤波器指标:wp=0.2πrad,ws=0.3πrad,αp=1dB,αs=40dB);

2. 用双线性变换法频率转换公式,将DF技术指标转换为巴特沃斯AF的技术指标;

3. 调用MATLAB函数buttord和butter,设计该模拟滤波器;

4. 根据1所确定的技术指标,调用MATLAB函数buttord和butter,直接设计数字滤波器,观察设计结果与上面用双线性变换法的设计结果是否相同。

5. 滤波仿真:调用MATLAB工具箱函数filter对下面给出的带加性干扰的信号x(n)进行滤波,观察滤波效果(绘制滤波前后信号的时域和频域波形)。a. 滤波仿真:调用MATLAB工具箱函数filter对下面给出的带加性干扰的信号进行滤波,观察滤波效果(绘制滤波前后信号的时域和频域波形)。

  实验信号x(n)的128个样值:

   -0.0289    0.3943    0.9965    1.1266    0.9995    1.0891    1.2262    1.0699

    0.8990    0.7685    0.7844    0.9471    1.4317    1.6765    1.7629   -3.2903

    3.4122    4.5403   -2.1667    -2.0584    4.6694    2.0368   -0.4864    0.1427

    2.5652   -1.8980   -0.0527    -1.4730    2.7884   -6.4092    2.9084   -5.1428

    1.5929    0.0014   -0.6010    -4.3059   -0.4518    1.9959   -3.3526    0.5745

   -3.5487    0.5913   -0.2472    -1.5479   -2.4422    2.5066   -4.2421   -2.3588

    3.8869   -3.9855    0.9583    -1.2164    3.7050    1.2411   -1.7249    0.9964

    3.9695    1.3400   -3.5513    5.8552   -2.8092    2.6877    0.4444    3.5641

   -2.3496    3.6065   -1.7598    1.4699    3.7201   -1.1626    2.5171   -1.8247

    2.5076    2.5423   -0.1164    0.3413   -2.8366   -0.5732    1.1705   -1.1689

   -1.8778   -2.8823    0.1312    0.1701   -0.3147   -3.2178   -0.2897   -0.2661

   -2.6156    2.5836   -6.1123    0.3611    3.8320   -5.1405   -2.4962    1.3115

    0.3410   -5.0227    2.5150   -2.2485    0.8696   -0.4257   -0.2326    5.2529

   -1.1389   -1.0966    4.2358   -0.8846    1.3454    2.1462   -0.2605    1.8075

   -1.5669    3.4374    2.4737   -0.9521    1.0952    0.2180    0.3156    4.8910

   -1.5372    1.4585   -1.0904    5.1795   -1.9366   -1.0818    1.2667   -0.4268

 

MATLAB函数简介

● filter: 一维数字滤波器直接Ⅱ型实现函数。

yn=filter(B,A,xn) 按直接Ⅱ型实现结构对输入信号向量xn滤波,返回滤波器输出信号向量yn,调用参数B和A分别为滤波器系统函数分子和分母多项式系数向量。其实质是求解差分方程

A(1)y(n) = B(1)x(n) + B(2)x(n-1) + ... + B(M+1)x(n-M) -A(2)y(n-1) - ... -A(N+1)y(n-N)

如果A(1)不等于1时,则对系数关于A(1)归一化后计算输出信号y(n)。当A=1时对应FIRDF的直接Ⅱ型实现

● freqz: 求数字滤波器频率响应。

[H,W]=freqz(B,A,N) 计算出数字滤波器频率向应函数在[0,π]上的N点等间隔采样复数向量H和N个采样频率点向量W。调用参数B和A分别为滤波器系统函数分子和分母多项式系数向量。

[H,W] = freqz(B,A,N,'whole') 计算出数字滤波器频率向应函数在[0,2π]上的N点等间隔采样复数向量H和N个采样频率点向量W。

freqs:求模拟滤波器频率响应.请用help 命令查其功能和调用各式。

%用双线性变换法设计DF
T=1;Fs=1/T;%T=1s为方便起见
wpz=0.2;wsz=0.3;
%实验步骤c.d中涉及的程序
wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=40;%化为模拟滤波器的指标
[N,wc]=buttord(wp,ws,rp,rs,'s');%设计该滤波器的过渡模拟滤波器
[B,A]=butter(N,wc,'s');
[Bz,Az]=bilinear(B,A,Fs);%用双线性变换法转换成数字滤波器
%实验步骤g中涉及的程序
[Nd,wdc]=buttord(wpz,wsz,rp,rs);%调用buttord和butter直接设计数字滤波器
[Bdz,Adz]=butter(Nd,wdc);
[Hz,wz]=freqz(Bz,Az);
[Hdz,wdz]=freqz(Bdz,Adz);
figure(1)
subplot(211);
plot(wz/pi,10*log(abs(Hz)));
xlabel('\omega/\pi'),ylabel('Magnitude/dB');grid on
title('双线性变换法的幅频响应')
subplot(212);
plot(wdz/pi,10*log(abs(Hdz)));
xlabel('\omega/\pi'),ylabel('Magnitude/dB');grid on
title('直接设计数字滤波器的相频响应');grid on
Bdz,Adz%调用buttord和butter直接设计数字滤波器的系统函数的分子分母多项式系数向量B和A
Bz,Az%双线性变换法数字滤波器系统函数的分子分母多项式系数向量B和A
N  %输出阶数,与自己在c算出来的进行比较,说明是正确的


%%%%%%%g步骤结论%%%%%%
%二者所设计的滤波器的幅频响应基本一致,但是直接设计的滤波器的过渡带的斜率的绝对值比
%双线性的小,即双线性的幅频响应较好,而且在w=pi(高频)处没有抖动的波纹

%%%%%让信号通过滤波器

xn=[-0.0289    0.3943    0.9965    1.1266    0.9995    1.0891    1.2262    1.0699...
    0.8990    0.7685    0.7844    0.9471    1.4317    1.6765    1.7629   -3.2903...
    3.4122    4.5403   -2.1667    -2.0584    4.6694    2.0368   -0.4864    0.1427...
    2.5652   -1.8980   -0.0527    -1.4730    2.7884   -6.4092    2.9084   -5.1428...
    1.5929    0.0014   -0.6010    -4.3059   -0.4518    1.9959   -3.3526    0.5745...
    -3.5487    0.5913   -0.2472    -1.5479   -2.4422    2.5066   -4.2421   -2.3588...
    3.8869   -3.9855    0.9583    -1.2164    3.7050    1.2411   -1.7249    0.9964...
    3.9695    1.3400   -3.5513    5.8552   -2.8092    2.6877    0.4444    3.5641...
    -2.3496    3.6065   -1.7598    1.4699    3.7201   -1.1626    2.5171   -1.8247...
    2.5076    2.5423   -0.1164    0.3413   -2.8366   -0.5732    1.1705   -1.1689...
    -1.8778   -2.8823    0.1312    0.1701   -0.3147   -3.2178   -0.2897   -0.2661...
    -2.6156    2.5836   -6.1123    0.3611    3.8320   -5.1405   -2.4962    1.3115...
    0.3410   -5.0227    2.5150   -2.2485    0.8696   -0.4257   -0.2326    5.2529...
    -1.1389   -1.0966    4.2358   -0.8846    1.3454    2.1462   -0.2605    1.8075...
    -1.5669    3.4374    2.4737   -0.9521    1.0952    0.2180    0.3156    4.8910...
    -1.5372    1.4585   -1.0904    5.1795   -1.9366   -1.0818    1.2667   -0.4268];

sigLength=length(xn);%信号长度
Y = fft(xn,sigLength);%做FFT变化
Pyy = Y.* conj(Y) / sigLength; %求出幅值
halflength=floor(sigLength/2); 
f=Fs*(0:halflength)/sigLength; 

y=filter(Bz,Az,xn);
k=fft(y);                      %滤波后的波形做离散傅里叶变换
w=2*[0:length(k)-1]/length(k);

figure(2);
subplot(211)
plot(f,Pyy(1:halflength+1))
xlabel('Frequency(Hz)'),ylabel('幅度')%画未经过滤波的频域波形
grid on

subplot(212)
plot(w,abs(k));
xlabel('w/pi')
ylabel('幅度')
title('IIR滤波器滤波后频谱');  %解调滤波后的频谱
grid on

xn0=ifft(y);%逆变换,恢复时域信号
xn0=abs(xn0)

figure(3)
subplot(211);
stem(xn0)   %plot(xn0);
xlabel('Time(s)');
title('经过滤波器信号时域波形');
grid on
subplot(212);
stem(xn)    %plot(xn0);
xlabel('Time(s)');
title('原始信号时域波形');
grid on
%%%%%%%h步骤的结论%%%%%%%%%
%无论是从时域还是频域都可以看出大于0..3的频率部分应经基本被滤除,效果较为理想

 

Category: DSP | Tags: matlab IIR | Read Count: 2056

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com