width=2*pi*1200/fs-2*pi*1000/fs; n=ceil(12.8*pi/width)+1; wn=(ws+wp)/2; beta=10.056;
win=kaiser(n,beta); b=fir1(n-1,wn/pi,win); [H,m]=freqz(b,[1],200); figure(1)
plot(m,abs(H));
title('FIR数字带通滤波器幅度响应') figure(2)
subplot(2,1,1) plot(t,y);
title('声音信号时域波形') c=fftfilt(b,y); subplot(2,1,2) plot(t,c);
title('过滤后声音信号时域波形') y1=fft(y); c1=fft(c); figure(3)
subplot(2,1,1);
plot(f,abs(y1(1:N/2))); title('声音信号幅度响应') subplot(2,1,2);
plot(f,abs(c1(1:N/2)));
title('过滤后声音信号幅度响应') sound(c,fs)
(4)IIR低通
function[]=iirlowpass();
[x,fs,bits]=wavread('2.wav');%播放原始信号 N=length(x);%返回采样点数 t=(1:N)/fs;
df=fs/N;%采样间隔 n1=1:N/2;
f=(n1-1)*df;%频带宽度
wp=2000/(fs/2);ws=2400/(fs/2); [n,wn]=buttord(wp,ws,1,10); [b,a]=butter(n,wn); [H,m]=freqz(b,a); figure(1);
plot(m*fs/(2*pi),abs(H));
title('IIR低通滤波器的频率响应'); xlabel('频率W(rad)');
ylabel('幅值'); grid; figure(2);
y3=fft(x);%原信号进行快速傅立叶变换
ys=filter(b,a,x);%信号送入滤波器滤波,ys为输出 yz=fft(ys); subplot(2,1,1);
plot(f,abs(y3(1:N/2))); title('原信号的频谱图'); xlabel('频率'); ylabel('幅值'); subplot(2,1,2);
plot(f,abs(yz(1:N/2)));
title('IIR低通滤波器滤波后的信号频谱图'); xlabel('频率');
ylabel('幅值');%声音信号时域波形 figure(3);
subplot(2,1,1); plot(t,x);
title('声音信号时域波形'); subplot(2,1,2); plot(t,ys);
title('滤波后声音信号时域波形'); sound(ys,fs); (5)IIR高通
function[]=iirhighpass(); [x,fs,bits]=wavread('2.wav'); N=length(x);%返回采样点数 t=(1:N)/fs;
df=fs/N;%采样间隔 n1=1:N/2;
f=(n1-1)*df;%频带宽度
wp=5000/(fs/2);ws=4700/(fs/2); [n,wn]=buttord(wp,ws,1,10); [b,a]=butter(n,wn,'high'); [H,m]=freqz(b,a); figure(1);
plot(m*fs/(2*pi),abs(H));
title('巴特沃斯滤波器的频率响应'); xlabel('频率W(rad)'); ylabel('幅值'); grid; figure(3);
ys=filter(b,a,x);%信号送入滤波器滤波,ys为输出
ya=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 y3=fft(x);%原信号进行快速傅立叶变换 subplot(2,1,1);
plot(f,abs(y3(1:N/2))); title('原信号的频谱图'); xlabel('频率'); ylabel('幅值'); subplot(2,1,2);
plot(f,abs(ya(1:N/2)));
title('IIR高通滤波器滤波后的信号频谱图'); xlabel('频率'); ylabel('幅值'); figure(2)
subplot(2,1,1); plot(t,x);
title('原信号时域图'); subplot(2,1,2); plot(t,ys);
title('滤波后信号时域图'); sound(ys,fs); (6)IIR带通
function[]=iirbandpass(); [x,fs,bits]=wavread('2.wav'); N=length(x);%返回采样点数 t=(1:N)/fs;
df=fs/N;%采样间隔 n1=1:N/2;
f=(n1-1)*df;%频带宽度
wp1=2400/(fs/2);ws1=2000/(fs/2); wp2=6000/(fs/2);ws2=6400/(fs/2);
[n,wn]=buttord([wp1 wp2],[ws1 ws2],2,10); [b,a]=butter(n,wn); [H,m]=freqz(b,a); figure(1);
plot(m*fs/(2*pi),abs(H));
title('巴特沃斯滤波器的频率响应'); xlabel('频率W(rad)'); ylabel('幅值'); grid; figure(2);
ys=filter(b,a,x);%信号送入滤波器滤波,ys为输出 ya=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 y3=fft(x);%原信号进行快速傅立叶变换 subplot(2,1,1);
plot(f,abs(y3(1:N/2))); title('原信号的频谱图'); xlabel('频率'); ylabel('幅值'); subplot(2,1,2);
plot(f,abs(ya(1:N/2)));
title('IIR带通滤波器滤波后的信号频谱图'); xlabel('频率'); ylabel('幅值'); figure(3);
subplot(2,1,1); plot(t,x);
title('原信号时域'); subplot(2,1,2); plot(t,ys);
title('滤波后信号时域'); sound(ys,fs);
实验八、有限冲击响应滤波器(FIR)算法实验
一、实验目的
1、掌握用窗函数法设计FIR数字滤波器的原理和方法。 2、熟悉线性相位FIR数字滤波器特性。 3、了解各种窗函数对滤波特性的影响。
二、实验设备
计算机、ZY13DSP12BC2实验箱。
三、实验原理
1、有限冲击响应数字滤波器的基础理论。
FIR数字滤波器是一种非递归系统,其冲激响应h(n)是有限长序列,其差分方程表达式为:
y(n)??h(i)x(n?i)
i?0N?1N为FIR滤波器的阶数。
在数字信号处理应用中往往需要设计线性相位的滤波器,FIR滤波器在保证幅度特性满足技术要求的同时,很容易做到严格的线性相位特性。为了使滤波器满足线性相位条件,要求其单位脉冲响应h(n)为实序列,且满足偶对称或奇对称条件,即h(n)=h(N-1-n)或h(n)=-h(N-1-n)。这样,当N为偶数时,偶对称线性相位FIR滤波器的差分方程表达式为
N/2?1 y(n)??h(i)(x(n?i)?x(N?1?n?i))
i?0由上可见,FIR滤波器不断地对输入样本x(n)延时后,再做乘法累加算法,将滤波器结果y(n)输出,因此,FIR实际上是一种乘法累加运算。而对于线性相位FIR而言,利用线性相位FIR滤波器系数的对称特性,可以采用结构精简的FIR结构将乘法器数目减少一半。
2、 模拟滤波器原理(巴特沃斯滤波器、切比雪夫滤波器)
为了用物理可实现的系统逼近理想滤波器的特性,通常对理想特性作如下修改:
(1)允许滤波器的幅频特性在通带和阻带有一定的衰减范围,幅频特性在这一范围内允许有起伏。
(2)在通带与阻带之间允许有一定的过渡带。
工程中常用的逼近方式有巴特沃思(Butterworth)逼近、切比雪夫(Chebyshev)逼近和椭圆函数逼近。相应设计的滤波器分别为巴特沃思滤波器、切比雪夫滤波器和椭圆函数滤波器。
巴特沃思滤波器的模平方函数由下式描述: |HB(?)|2?1 n为阶数;?c为滤波器截止频率 ?2n1?()?c