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