数字滤波器是一个离散的系统（按预定的算法，将输入离散时间信号转换为所要求的输出离散时间信号的特定功能装置）。与模拟滤波器相比，数字滤波器具有很多突出的优点，例如它可以满足滤波器对幅度和相位特性的严格要求，可以避免模拟滤波器所无法克服的电压漂移和噪声问题。数字滤波器分为有限冲激响应数字滤波器，即FIR数字滤波器和无限冲激响应，即IIR数字滤波器。在数字信号处理中，利用数字滤波器可改变信号中所含频率分量的相对比例或滤除某些频率分量，使其达到所需要的效果．其中数字FIR滤波器由于具有精确的线性相位，且系统稳定，所以广泛应用于通 信、数字图象处理、语音信号处理、自适应处理、雷达／声纳系统等方面．目前，FIR滤波器的设计主要有窗函数设计法和频率采样设计法．但是这两种方法都不易精确控制通带边界频率与阻带边界频率，所以在实际应用中具有一定的局限性而以最大误差最小化准则支持的切比雪夫逼近法是一种优异的设计方法，易于精确控制wp与ws。
等波纹滤波器的最优化设计方法主要有2种，第1种是离散最小二乘法．它的思路是在给定的一些离散点上，使实际的幅频特性和理想幅频特性之间的误差的平方和为最小．第2种是最优化等波纹设计法，也称为雷米兹法或切比雪夫逼近法．该类型滤波器幅频特性在通带和阻带上的误差峰值是均匀分布的，其误差具有等波纹特性，因而把波纹的幅度控制到最小，或在同等指标下减小它的阶次．第1种方法是连续最小的平方法的推广，容易理解，但它的指标与滤波器没有直接关联，误差平方小的滤波器不能保证没有窄而大的波纹出现，像吉布斯效应那样．第2种方法直接控制通带波动和阻带衰减，最具针对性，是滤波器的最优化设计方法．因此，采用MATLAB 信号处理工具箱提供的函数，运用最优化等波纹设计法实现数字FIR滤波器的设计和仿真．
完整的最优化等波纹滤波器设计，除了切比雪夫等波纹逼近公式外，还要考虑：
（1）滤波器采样点数n 的确定；
（2）)极值数目的确定．最优化等波纹滤波器的误差函数在给定的频率上有(L+2)或(L+3)个极值，
L为多项式的阶数．对于某些Wp，Ws的组合可能得到有(L+3)个极值的滤波器．
(3) 建立频率修正的算法．在程序中自动进行反复的迭代修正，直到达到要求的精度为止．交替定理能保证切比雪夫逼近问题的解存在并且惟一，但它并没有说明如何得到该解，既不知n(或L)，也不知极值的频率Wi 和波纹系数
, Parks和McClellan提供了利用Remez交换算法导出的迭代算法。
决定切比雪夫等波纹逼近低通滤波器系数的参数主要有：滤波器长度M，通带和阻带截止频率wp、ws，相应频带的幅度m，权系数w。其中权系数w由通带和阻带波动Ap、Ar决定。使用权系数w，是考虑在设计滤波器时对通带和阻带常要求不同的逼近精度，故乘以不同的权系数，以统一使用最小化最大误差。
在优化设计的Matlab实现中，程序中经常使用remez函数，这种函数的使用方法为： b=remez(n,f,a,w,’ftype’)
①n为待设计滤波器的阶数；f是一个向量，它是一个0到1的正数
②a是一个向量，指定频率段的幅度值；w对应于各个频段的加权值
③函数的返回值b是设计出的滤波器的系数组成的一个长度为n+1的向量
(1) 利用Remez函数设计等波纹低通滤波器
设计要求：①通带截频0.5，阻带截频0.6 , 采样频率2000Hz
②阻带衰减大于等于40dB，通带波纹0.1710和阻带波纹0.01
程序参见附录二中的3-（1）利用Remez函数设计等波纹低通滤波器
1.窗函数法的Matlab实现
(1) 利用fir1函数及hamming窗设计低通滤波器
f1=100;f2=200;%待滤波正弦信号频率
fs=2000;%采样频率
m=(0.3*f1)/(fs/2);%定义过度带宽
M=round(8/m);%定义窗函数的长度
N=M-1;%定义滤波器的阶数
b=fir1(N,0.5*f2/(fs/2));%使用fir1函数设计滤波器
%输入的参数分别是滤波器的阶数和截止频率
figure(1)
[h,f]=freqz(b,1,512);%滤波器的幅频特性图
%[H,W]=freqz(B,A,N)当N是一个整数时函数返回N点的频率向量和幅频响应向量
plot(f*fs/(2*pi),20*log10(abs(h)))%参数分别是频率与幅值
xlabel('频率/赫兹');ylabel('增益/分贝');title('滤波器的增益响应');
figure(2)
subplot(211)
t=0:1/fs:0.5;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%滤波前信号
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);%将信号变换到频域
AFs=abs(Fs);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s);%使用filter函数对信号进行滤波
%参数分别为滤波器系统函数的分子和分母多项式系数向量和待滤波信号输入
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
axis([0.2 0.5 -2 2]);%限定图像坐标范围
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
（2）利用fir1函数及Kaiser窗设计带通滤波器
Rs=0.01;fs=8000;%采样频率
fcuts=[1000 1300 2210 2410];a=[0 1 0];
dev=Rs*ones(1,length(a));
[M,Wc,beta,ftype]=kaiserord(fcuts,a,dev,fs);
%M为能够满足设计要求的滤波器的最小阶数，Wc为滤波器的截止频率点
%第一个元素f为待设计滤波器的过渡带的起始点和结束点
%第二个元素a指定第一个元素频率段的理想幅度值
%第三个元素dev中的元素为各通带和阻带内允许的幅度最大误差
M=mod(M,2)+M;
window=Kaiser(M+1,beta);b=fir1(M,Wc,ftype,window);
%输入的第一个参数是滤波器的阶数
%第二个参数是滤波器的截止频率
%第三个参数是滤波器的类型，stop为带阻滤波器
%第四个参数是采用的窗函数
[h,f]=freqz(b,1,512);%滤波器的幅频特性图
%[H,W]=freqz(B,A,N)当N是一个整数时函数返回N点频率向量和幅频响应向
figure(1)
plot(f*fs/(2*pi),20*log10(abs(h)))%参数分别是频率与幅值
xlabel('频率/赫兹');ylabel('增益/分贝');title('滤波器的增益响应');
f1=500;f2=1500;f3=2000;f4=3000;%待滤波正弦信号频率
t=(0:200)/fs;%定义时间的步长
s=sin(2*f1*pi*t)+sin(2*f2*pi*t)+sin(2*f3*pi*t)+sin(2*f4*pi*t);
sf=filter(b,1,s);%使用filter函数对信号进行滤波
figure(2)
subplot(211)
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255);
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
axis([0.005 0.025 -4 4]);
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
（3）利用fir1函数及Kaiser窗设计多通带滤波器
Rs=0.01;fs=200;%采样频率
fcuts=[10 20 40 50 60 70 80 90];a=[0,1,0,1,0];
dev=Rs*ones(1,length(a));
[M,Wc,beta,ftype]=kaiserord(fcuts,a,dev,fs);
M=mod(M,2)+M;window=Kaiser(M+1,beta);b=fir1(M,Wc,ftype,window);
[h,f]=freqz(b,1,512);%滤波器的幅频特性图
figure(1)
plot(f*fs/(2*pi),20*log10(abs(h)))%参数分别是频率与幅值
xlabel('频率/赫兹');ylabel('增益/分贝');title('滤波器的增益响应');
f1=5;f2=20;f3=30;f4=55;f5=75;f6=95;%待滤波正弦信号频率
t=(0:200)/fs;%定义时间的步长
s1=sin(2*f1*pi*t)+sin(2*f2*pi*t)+sin(2*f3*pi*t);
s=s1+sin(2*f4*pi*t)+sin(2*f5*pi*t)+sin(2*f6*pi*t);%滤波前信号
sf=filter(b,1,s);%使用filter函数对信号进行滤波
figure(2)
subplot(211)
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255);
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
subplot(212)
Fsf=fft(sf,512); AFsf=abs(Fsf);%滤波后的信号频域图及信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图')
2.频率取样法的Matlab实现
(1) 利用频率取样法设计低通滤波器
M=63; Wp=0.5*pi;%所需频率采样点个数及通带截止频率
m=0:(M+1)/2; Wm=2*pi*m./(M+1);%通频带上的采样点及阻带截止频率
mtr=floor(Wp*(M+1)/(2*pi))+2;%向负方向入floor(3.5)=3;floor(-3.2)=-4
Ad=[Wm<=Wp];Ad(mtr)=0.38;
Hd=Ad.*exp(-j*0.5*M*Wm);%构造频域采样向量H(k)
Hd=[Hd conj(fliplr(Hd(2:(M+1)/2)))];
%fliplr函数实现矩阵的左右翻转conj是求复数的共轭
h=real(ifft(Hd));%h（n）=IDFT[H(k)]
w=linspace(0,pi,1000);%用于产生0,pi之间的1000点行矢量
H=freqz(h,[1],w);%滤波器的幅频特性图
figure(1)
plot(w/pi,20*log10(abs(H)));%参数分别是归一化频率与幅值
xlabel('归一化角频率');ylabel('增益/分贝');title('滤波器的增益响应');
axis([0 1 -50 0.5]);
f1=100;f2=300;f3=700;fs=2000;%待滤波正弦信号频率及采样频率
figure(2)
subplot(211)
t=0:1/fs:0.25;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%滤波前信号
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512); AFs=abs(Fs);%将信号变换到频域及信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(h,1,s);%使用filter函数对信号进行滤波
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
axis([0.2 0.25 -2 2]);%限定图像坐标范围
subplot(212)
Fsf=fft(sf,512); AFsf=abs(Fsf);%滤波后的信号频域图及信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
（2）利用频率取样法设计高通滤波器
M=32;%所需频率采样点个数
Wp=0.6*pi;%通带截止频率
m=0:M/2;%阻频带上的采样点
Wm=2*pi*m./(M+1);%阻带截止频率
mtr=ceil(Wp*(M+1)/(2*pi));%向正方向舍入ceil(3.5)=4;ceil(-3.2)=-3;
Ad=[Wm>=Wp];
Ad(mtr)=0.28;
Hd=Ad.*exp(-j*0.5*M*Wm);%构造频域采样向量H(k)
Hd=[Hd conj(fliplr(Hd(2:M/2+1)))];
%fliplr函数实现矩阵的左右翻转conj是求复数的共轭
h=real(ifft(Hd));%h（n）=IDFT[H(k)]
w=linspace(0,pi,1000);%用于产生0,pi之间的1000点行矢量
H=freqz(h,[1],w);%滤波器的幅频特性图
figure(1)
plot(w/pi,20*log10(abs(H)));%参数分别是归一化频率与幅值
xlabel('归一化频率');ylabel('增益/分贝');title('滤波器的增益响应');
axis([0 1 -50 0]);
f1=200;f2=700;f3=800%待滤波正弦信号频率
fs=2000;%采样频率
figure(2)
subplot(211)
t=0:1/fs:0.25;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%滤波前信号
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);%将信号变换到频域
AFs=abs(Fs);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(h,1,s);%使用filter函数对信号进行滤波
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
axis([0.2 0.25 -2 2]);%限定图像坐标范围
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
3.优化设计的Matlab实现
(1) 利用Remez函数设计等波纹低通滤波器
fs=2000;%设定采样频率
rp=3;%通带波纹
rs=40;%阻带波纹
f=[500 600];%截止频率
a=[1 0];%期望幅度
dev=[(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[n,fo,ao,w]=remezord(f,a,dev,fs);
%函数remezord返回参数n表示滤波器的阶数
%FIR滤波器有B个频带时，f,a,dev分别为2B-2，B,B个元素的向量
b=remez(n,fo,ao,w);%函数remez的返回值为n阶FIR滤波器的系数
%fo，ao是2B个元素的向量，分别表示B个频带的2B边界频率及幅度值
%w是B个元素的向量，表示各频带的加权值
figure(1)
freqz(b,1,1024,fs);%滤波器的特性图
f1=400;f2=700;%待滤波正弦信号频率
t=0:1/fs:0.1;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%滤波前信号
figure(2)
subplot(211)
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);%将信号变换到频域
AFs=abs(Fs);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s);%使用filter函数对信号进行滤波
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
(2) 利用Remez函数设计等波纹带通滤波器
fs=2000;%设定采样频率
rp=3;%通带波纹
rs=40;%阻带波纹
f=[200 300 600 700];%截止频率
a=[0 1 0];%期望幅度
dev=[10^(-rs/20) (10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[n,fo,ao,w]=remezord(f,a,dev,fs);
b=remez(n,fo,ao,w);
figure(1)
freqz(b,1,1024,fs);%滤波器的特性图
f1=100;f2=400;f3=500;f4=800;%待滤波正弦信号频率
t=0:1/fs:0.1;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t)+sin(2*pi*f4*t);
%滤波前信号
figure(2)
subplot(211)
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);%将信号变换到频域
AFs=abs(Fs);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s);%使用filter函数对信号进行滤波
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
(3) 利用Remez函数设计等波纹带阻滤波器
fs=2000;%设定采样频率
rp=3;%通带波纹
rs=40;%阻带波纹
f=[200 300 600 700];%截止频率
a=[1 0 1];%期望幅度
dev=[10^(-rs/20) (10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[n,fo,ao,w]=remezord(f,a,dev,fs);
b=remez(n,fo,ao,w);
figure(1)
freqz(b,1,1024,fs);%滤波器的特性图
f1=100;f2=400;f3=500;f4=800;%待滤波正弦信号频率
t=0:1/fs:0.1;%定义时间范围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t)+sin(2*pi*f4*t);%滤波前信号
figure(2)
subplot(211)
plot(t,s);%滤波前的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512);%将信号变换到频域
AFs=abs(Fs);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFs(1:256));%滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s);%使用filter函数对信号进行滤波
subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图');
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
4.利用滤波器处理加有噪声的音频波形
- 利用窗函数法设计的低通滤波器处理加有噪声的音频波形
f3=7000;%所加噪声正弦函数的频率
[Y,fs,bits]=wavread('D:\yinpin\nihao01.wav');
%利用wavread产生音频的函数及采样频率
L=length(Y);t=0:1/fs:(L-1)/fs;%定义时间的范围及步长
y=0.005*sin(2*pi*f3*t); n1=floor(L/2);%所加噪声
f1=(0:n1)*fs/L;Y=Y(:,1);sound(Y,fs);%输出加噪前音频
Y1=y+Y';%给音频加噪声
FY1=abs(fft(Y1,L));FY=abs(fft(Y,L));sound(Y1,fs);%输出加噪后的音频
figure(1)
subplot(211)
plot(t(1:1000),Y(1:1000)); grid on;%加噪前音频的时域图
xlabel('时间(t)');ylabel('幅度(Y)');
title('加噪前录音波形的时域图');
subplot(212)
plot(f1,FY(1:n1+1)); grid on;%加噪前音频的频域图
xlabel('频率(f)');ylabel('幅度(FY)');
title('加噪前录音波形的频域图');
figure(2)
subplot(211)
plot(t(1:1000),Y1(1:1000)); grid on;%加噪后音频的时域图
xlabel('时间(t)');ylabel('幅度(Y1)');
title('加噪声后录音波形的时域图');
subplot(212)
plot(f1,FY1(1:n1+1)); grid on;%加噪后音频的频域图
xlabel('频率(f)');ylabel('幅度(FY1)');
title('加噪声后录音波形的频域图');
m=0.03; M=round(8/m);
N=M-1;%定义滤波器的阶数
b=fir1(N,0.6);
figure(3)
[h,f]=freqz(b,1,512);%滤波器的幅频特性图
plot(f*fs/(2*pi),20*log10(abs(h)))%参数分别是频率与幅值
xlabel('频率/赫兹');ylabel('增益/分贝');title('滤波器的增益响应');
figure(4)
sf=filter(b,1,Y1);%使用filter函数对信号进行滤波
Fsf=abs(fft(sf,L));
subplot(211)
plot(t(1:1000),sf(1:1000)); grid on;%滤波后音频的时域图
xlabel('时间(t)');ylabel('幅度(sf)');
title('滤波后录音波形的时域图');
axis([0.01 0.05 -0.002 0.002])
subplot(212)
plot(f1,Fsf(1:n1+1)); grid on;%滤波后音频的频域图
xlabel('频率(f)');ylabel('幅度(Fsf)');
title('滤波后录音波形的频域图');sound(sf,fs);
- 利用优化设计的低通滤波器处理加有噪声的音频波形
f3=7000;%所加噪声正弦函数的频率
[Y,fs,bits]=wavread('D:\yinpin\nihao01.wav');
%利用wavread产生音频的函数及采样频率
L=length(Y);t=0:1/fs:(L-1)/fs;%定义时间的范围及步长
y=0.005*sin(2*pi*f3*t);%所加噪声
n1=floor(L/2);f1=(0:n1)*fs/L;Y=Y(:,1);sound(Y,fs);%输出加噪前音频
Y1=y+Y';%给音频加噪声
FY1=abs(fft(Y1,L));FY=abs(fft(Y,L));sound(Y1,fs);%输出加噪后的音频
figure(1)
subplot(211)
plot(t(1:1000),Y(1:1000)); grid on;%加噪前音频的时域图
xlabel('时间(t)');ylabel('幅度(Y)');
title('加噪前录音波形的时域图');
subplot(212)
plot(f1,FY(1:n1+1)); grid on;%加噪前音频的频域图
xlabel('频率(f)');ylabel('幅度(FY)');
title('加噪前录音波形的频域图');
figure(2)
subplot(211)
plot(t(1:1000),Y1(1:1000)); grid on;%加噪后音频的时域图
xlabel('时间(t)');ylabel('幅度(Y1)');
title('加噪声后录音波形的时域图');
subplot(212)
plot(f1,FY1(1:n1+1)); grid on;%加噪后音频的频域图
xlabel('频率(f)');ylabel('幅度(FY1)');
title('加噪声后录音波形的频域图');
rp=3; rs=40; %通带波纹及阻带波纹
f=[6000 7000];%截止频率
a=[1 0];%期望幅度
dev=[(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)];
[n,fo,ao,w]=remezord(f,a,dev,fs); b=remez(n,fo,ao,w);
figure(3)
freqz(b,1,1024,fs);%滤波器的特性图
figure(4)
sf=filter(b,1,Y1);%使用filter函数对信号进行滤波
Fsf=abs(fft(sf,L));
subplot(211)
plot(t(1:1000),sf(1:1000));%滤波后音频的时域图
grid on;
xlabel('时间(t)');ylabel('幅度(sf)');
title('滤波后录音波形的时域图');
subplot(212)
plot(f1,Fsf(1:n1+1));%滤波后音频的频域图
grid on;
xlabel('频率(f)');ylabel('幅度(Fsf)');
title('滤波后录音波形的频域图');
sound(sf,fs);
图4-1 等波纹低通滤波器的增益响应
从参考程序及图4-1可以得到所设计出滤波器的参数如下：
①滤波器的采样频率为2000Hz，滤波器的阶数为22
②滤波器的通带截频0.5，阻带截频0.6，过渡带宽均为0.1
③阻带衰减为40dB，通带波纹为0.1710，阻带波纹为0.01
对比设计要求与所设计出滤波器的参数可知，其各项参数均满足设计指标，所设计出的滤波器即为设计所要求的滤波器。
图4-2 信号滤波前的时域图和频域图