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