设计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的频率部分应经基本被滤除，效果较为理想