基于matlab 的数字滤波器
 


clear
clc
fs = 22050;
wp=0.1*pi;%通带截止频率
ws=0.4*pi;%阻带截止频率
Rp=3;%通带衰减率
Rs=75;%阻带衰减率
Fs=22050;
Ts=1/Fs;
wp1=2/Ts*tan(wp/2);                 %将模拟指标转换成数字指标
ws1=2/Ts*tan(ws/2); 
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');  %选择滤波器的最小阶数，求出3db的截止频率
[Z,P,K]=buttap(N);                  %创建butterworth模拟滤波器，零点，极点，增益
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);   
[bd,ad]=bilinear(b,a,Fs);      %用双线性变换法实现模拟滤波器到数字滤波器的转换
[h,w]=freqz(bd,ad); 
figure(1)
subplot(111);
plot(w*fs/(2*pi),abs(h))
grid;
title('滤波器的性能分析'); 
figure(2)
[x,fs]=audioread('wslhjz.wav');
n=length(x);
f=fs*(0:fix(n/2-1))/n;%归一化
X=fft(x);
z=filter(bd,ad,x);
subplot(211);
plot(x);
title('原始信号的波形');
subplot(212);
plot(z);
title('滤波后信号的波形'); 
figure(3)
sound(z,fs);
subplot(211);
plot(f,abs(X(1:fix(n/2))));
title('原始信号的频谱');
xlabel('Hz');
Z=fft(z);
subplot(212);
plot(f,abs(Z(1:fix(n/2)))); 
title('滤波后的信号频谱');
xlabel('Hz');

在命令行中输入：
 
designfilt
 
然后就会弹出滤波器的种类，选定后会让你填参数。
 参数填好之后点确定，会在命令行窗口生成这个滤波器，复制到编辑器里直接用就好。
 
%读取音频文件%
[x,Fs]=audioread('C:\Users\16000\Desktop\testfile.wav');
%x是一个序列
%Fs是这个音频的采样率
%复制过来的滤波器%
High_Filt=designfilt('highpassfir', 'StopbandFrequency', 500, 'PassbandFrequency', 600, 'StopbandAttenuation', 60, 'PassbandRipple', 1, 'SampleRate', 48000);

%改成单声道%
x=x(:,1);
%滤波操作%
y_1 = filter(High_Filt,x);

sigLength=length(y_1); 
t=(0:length(y_1)-1)/Fs; %时间轴
subplot(2,1,1);
plot(t,y_1);xlabel('Time(s)');%在第一个窗口画波形
grid on;

 
使用designfilt方法的好处是自动生成，而且不用在filter命令里那里给滤波器填参数。
 
%部分重要参数说明%
Fs=Sample Rate;
Passband frequency是滤波器起始点
Stopband frequency是滤波器结束点
例如，Passband=400;Stopband=500的highpassfir就是说在400Hz到500Hz之间把波滤干净。
 
最后总结一下:
 
在命令行中输入designfilt;
设置好后将生成的代码复制进编辑器区域;
使用audioread(‘路径’)读入音频，拿到信号x和采样率Fs;
y=Filter(‘滤波器名称’，x)；
画图
 
一些其他看起来有点复杂的东西：
 
help filterDesigner		%滤波器涉及图形化界面
help butter	%butterworth滤波器

 
摘要
 
一、引言
 
二、正文
 
 1.设计要求
 
 2.设计步骤
 
 3.设计内容
 
 4.简易GUI设计
 
三、结论
 
四、收获与心得
 
五、附录
 
一、引言
 
随着Matlab仿真技术的推广，我们可以在计算机上对声音信号进行处理，甚至是模拟。通过计算机作图，采样，我们可以更加直观的了解语音信号的性质，通过matlab编程，调用相关的函数，我们可以非常方便的对信号进行运算和处理。
 


 
二、正文
 
2.1 设计要求
 
在有噪音的环境中录制语音，并设计滤波器去除噪声。
 
2.2 设计步骤
 
分析原始信号，画出原始信号频谱图及时频图，确定滤波器类型及相关指标 ；
按照类型及指标要求设计出滤波器，画出滤波器幅度和相位响应，分析该滤波器是否符合要求；
用所设计的滤波器对原始信号进行滤波处理，画出滤波后信号的频谱图及时频图；
对滤波前的信号进行分析比对，评估所设计滤波器性能。
 
2.3 设计内容
 
1.原始信号分析
 

 

  
 
 

 

  
 
 

 

  
 
 
分析信号的谱图可知，噪音在1650HZ和3300HZ附近的能量较高，而人声的能量基本位于1000HZ以下。因此，可以设计低通滤波器对信号进行去噪处理。
 
2.IIR滤波器设计
 
 用双线性变换法分别设计了巴特沃斯低通滤波器和椭圆低通滤波器和带阻滤波器：
 
①巴特沃斯滤波器
 
fp=800;fs=1300;rs=35;rp=0.5;
 
程序代码如下：
 
fp=800;fs=1300;rs=35;rp=0.5;Fs=44100;
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=buttord(wp,ws,rp,rs,'s');
 
[b,a]=butter(n,wn,'s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 


 

 

  
 
 

 

  
 
 

 

  
 
 


 


 
②椭圆低通滤波器
 
fp=1300;fs=1600;rs=60;rp=0.5;
 
程序代码如下：
 
fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=ellipord(wp,ws,rp,rs,'s');
 
[b,a]=ellip(n,rp,rs,wn,'s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 


 


 

 

  
 
 

 

  
 
 

 

  
 
 


 
③带阻滤波器
 
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6
 
fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6
 
程序代码如下：
 
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6;Fs=44100;
 
fp=[fp1,fp2];fs=[fs1,fs2];
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));
 
ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=buttord(wp,ws,rp,rs,'s');
 
[b,a]=butter(n,wn,'stop','s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 
fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6;Fs=44100;
 
fp1=[fp3,fp4];fs1=[fs3,fs4];
 
wp1=2*Fs*tan(2*pi*fp1/(2*Fs));
 
ws1=2*Fs*tan(2*pi*fs1/(2*Fs));
 
[n1,wn1]=buttord(wp1,ws1,rp,rs,'s');
 
[b1,a1]=butter(n1,wn1,'stop','s');
 
[num1,den1]=bilinear(b1,a1,Fs);
 
[h1,w1]=freqz(num1,den1,512,Fs);
 


 

 

  
 
 

 

  
 
 

 

  
 
 

 

  
 
 

 

  
 
 

 

  
 
 


 


 
FIR滤波器
 
①加hamming窗
 
n=100;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),Hamming(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 


 

 

  
 
 

 

  
 
 

 

  
 
 


 
②加hanning窗
 
n=;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),Hanning(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 


 

 

  
 
 

 

  
 
 

 

  
 
 


 


 
③加blackman窗
 
n=100;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),blackman(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 


 

 

  
 
 

 

  
 
 

 

  
 
 


 


 
滤波前后比对
 
 ①巴特沃斯低通滤波器滤波后
 

 

  
 
 

 

  
 
 

 

  
 
 


 
②椭圆低通滤波器滤波后
 


 

 

  
 
 

 

  
 
 

 

  
 
 
③带阻滤波器
 

 

  
 
 

 

  
 
 

 

  
 
 
④加hamming窗
 


 

 

  
 
 


 
⑤加hanning窗
 


 

 

  
 
 


 
⑥加blackman窗 
 


 

 

  
 
 


 
2.4简易GUI界面设计
 
 为了便于操作和演示，设计了如下的简易GUI界面。
 


 

 

  
 
 


 


 

 

  
 
 


 


 
结论
 
 由以上谱图分析可知，经过滤波器滤波后，信号中的高频杂音明显被抑制，而人声成分大部分被保留，起到了预期的滤波作用。
 
 对比所设计的两种滤波器，椭圆滤波器在过渡带相对较窄的情况下，能满足相对较高阻带衰减。
 
四、收获与心得
 
本次设计大概进行了五周的时间，语音信号处理的是目前比较流行且十分有趣的，在编程实现的过程中还是遇到了很多困难。在前期的准备工作中，查阅了大量资料，以完善我们的理论知识。我们为了完成本次设计，我们通过查阅相关书籍以及matlab中的帮助，选用不同的matlab函数，尝试不同的参数。经过接近一个礼拜的反复调试，最终基本的实现了设计任务。
 
虽然遇到了很多困难，但是我们在设计过程中都有收获很大。本次设计将信号与系统课上学习的知识用于实践，让我们对对语音信号处理更深入的了解，也让我们加深了对滤波器相关内容的理解，同时也使得我们的Matlab能力有了很大的提高。
 


 


 


 


 
参考文献
 
《应用matlab实现信号分析和处理》 科学出版社 
 
附录 
 
1 巴特沃斯低通滤波器
 
fp=800;fs=1300;rs=35;rp=0.5;Fs=44100;
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));
 
ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=buttord(wp,ws,rp,rs,'s');
 
[b,a]=butter(n,wn,'s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 
figure(1)
 
%subplot(3,1,1)
 
plot(w,abs(h));
 
xlabel('频率/Hz');ylabel('幅值');
 
title('巴特沃斯低通滤波器幅度特性');
 
axis([0,5000,0,1.2])
 
grid on;
 
figure(2)
 
%subplot(3,1,2)
 
plot(w,20*log10(abs(h)));
 
xlabel('频率/Hz');ylabel('幅值db');
 
title('巴特沃斯低通滤波器幅度特性db');
 
axis([0,5000,-90,10]);
 
grid on;
 
figure(3)
 
plot(w,180/pi*unwrap(angle(h)));
 
xlabel('频率/Hz');ylabel('相位');
 
title('巴特沃斯低通滤波器相位特性');
 
axis([0,5000,-1000,10])
 
grid on;
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);sound(x1,Fs,bits);
 
N1=length(x1);
 
Y1=fft(x1,N1);
 
f1=Fs*(0:N1-1)/N1;
 
t1=(0:N1-1)/Fs;
 
figure(4)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
y=filter(num,den,x1);
 
sound(y,Fs,bits);
 
N2=length(y);
 
Y2=fft(y,N2);
 
f2=Fs*(0:N2-1)/N2;
 
t2=(0:N2-1)/Fs;
 
figure(5)
 
plot(f2,abs(Y2))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 
2椭圆低通滤波器
 


 
fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));
 
ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=ellipord(wp,ws,rp,rs,'s');
 
[b,a]=ellip(n,rp,rs,wn,'s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 
figure(1)
 
plot(w,abs(h));
 
xlabel('频率/Hz');ylabel('幅值');
 
title('椭圆低通滤波器幅度特性');
 
axis([0,5000,0,1.2])
 
grid on;
 
figure(2)
 
plot(w,20*log10(abs(h)));
 
xlabel('频率/Hz');ylabel('幅值db');
 
title('椭圆低通滤波器幅度特性db');
 
axis([0,5000,-90,10]);
 
grid on;
 
figure(3)
 
plot(w,180/pi*unwrap(angle(h)));
 
xlabel('频率/Hz');ylabel('相位');
 
title('椭圆低通滤波器相位特性');
 
axis([0,5000,-1000,10])
 
grid on;
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);sound(x1,Fs,bits);
 
N1=length(x1);
 
Y1=fft(x1,N1); %对信号做N点FFT变换
 
f1=Fs*(0:N1-1)/N1;
 
t1=(0:N1-1)/Fs;
 
figure(4)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
y=filter(num,den,x1); 
 
sound(y,Fs,bits);
 
N2=length(y);
 
Y2=fft(y,N2); %对信号做N点FFT变换
 
f2=Fs*(0:N2-1)/N2;
 
t2=(0:N2-1)/Fs;
 
figure(5)
 
plot(f2,abs(Y2))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 


 
3.带阻滤波器
 
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6;Fs=44100;
 
fp=[fp1,fp2];fs=[fs1,fs2];
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));
 
ws=2*Fs*tan(2*pi*fs/(2*Fs));%wap=2*tan(wp/2)/Ts
 
[n,wn]=buttord(wp,ws,rp,rs,'s');
 
[b,a]=butter(n,wn,'stop','s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 
fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6;Fs=44100;
 
fp1=[fp3,fp4];fs1=[fs3,fs4];
 
wp1=2*Fs*tan(2*pi*fp1/(2*Fs));
 
ws1=2*Fs*tan(2*pi*fs1/(2*Fs));%wap=2*tan(wp/2)/Ts
 
[n1,wn1]=buttord(wp1,ws1,rp,rs,'s');
 
[b1,a1]=butter(n1,wn1,'stop','s');
 
[num1,den1]=bilinear(b1,a1,Fs);
 
[h1,w1]=freqz(num1,den1,512,Fs);
 
figure(1)
 
plot(w,abs(h));
 
xlabel('频率/Hz');ylabel('幅值');
 
title('巴特沃斯带阻滤波器幅度特性');
 
axis([0,5000,0,1.2])
 
grid on;
 
figure(2)
 
plot(w,20*log10(abs(h)));
 
xlabel('频率/Hz');ylabel('幅值db');
 
title('巴特沃斯带阻滤波器幅度特性db');
 
axis([0,5000,-90,10]);
 
grid on;
 
figure(3)
 
plot(w,180/pi*unwrap(angle(h)));
 
xlabel('频率/Hz');ylabel('相位');
 
title('巴特沃斯带阻滤波器相位特性');
 
axis([0,5000,-1000,10])
 
grid on;
 
figure(4)
 
plot(w1,abs(h1));
 
xlabel('频率/Hz');ylabel('幅值');
 
title('巴特沃斯带阻滤波器幅度特性');
 
axis([0,5000,0,1.2])
 
grid on;
 
figure(5)
 
plot(w1,20*log10(abs(h1)));
 
xlabel('频率/Hz');ylabel('幅值db');
 
title('巴特沃斯带阻滤波器幅度特性db');
 
axis([0,5000,-90,10]);
 
grid on;
 
figure(6)
 
plot(w1,180/pi*unwrap(angle(h1)));
 
xlabel('频率/Hz');ylabel('相位');
 
title('巴特沃斯带阻滤波器相位特性');
 
axis([0,5000,-1000,10])
 
grid on;
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);sound(x1,Fs,bits);
 
N1=length(x1);
 
Y1=fft(x1,N1); %对信号做N点FFT变换
 
f1=Fs*(0:N1-1)/N1;
 
t1=(0:N1-1)/Fs;
 
figure(7)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
y1=filter(num,den,x1);
 
y=filter(num1,den1,y1);
 
sound(y,Fs,bits);
 
N2=length(y);
 
Y2=fft(y,N2);
 
f2=Fs*(0:N2-1)/N2;
 
t2=(0:N2-1)/Fs;
 
figure(8)
 
plot(f2,abs(Y2))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 
加hamming窗
 
n=100;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),Hamming(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 
a=num2str(a);
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);
 
sound(x1,Fs,bits);
 
N1=length(x1);
 
Y1=fft(x1,N1);
 
f1=Fs*(0:N1-1)/N1;
 
y=fftfilt(b,x1);
 
sound(y,Fs,bits);
 
Y2=fft(y,N2);
 
f2=Fs*(0:N2-1)/N2;
 
figure(4)
 
subplot(2,1,1)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
subplot(2,1,2)
 
plot(f2,abs(Y2));
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 
加hanning窗
 
n=100;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),Hanning(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 
a=num2str(a);
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);
 
N1=length(x1);
 
Y1=fft(x1,N1);
 
f1=Fs*(0:N1-1)/N1;
 
y=fftfilt(b,x1);
 
sound(y,Fs,bits);
 
N2=length(y);
 
Y2=fft(y,N2);
 
f2=Fs*(0:N2-1)/N2;
 
figure(4)
 
subplot(2,1,1)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
subplot(2,1,2)
 
plot(f2,abs(Y2));
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 
加blackman窗
 
n=100;fp=1000;Fs=44100;
 
b=fir1(n,fp/(Fs/2),blackman(n+1));
 
[h,w]=freqz(b,1,512,Fs);
 
a=num2str(a);
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x1=s1(:,1);
 
sound(x1,Fs,bits);
 
N1=length(x1);
 
Y1=fft(x1,N1);
 
f1=Fs*(0:N1-1)/N1;
 
y=fftfilt(b,x1);
 
sound(y,Fs,bits);
 
N2=length(y);
 
f2=Fs*(0:N2-1)/N2;
 
figure(4)
 
subplot(2,1,1)
 
plot(f1,abs(Y1))
 
xlabel('频率/Hz');ylabel('幅度');
 
title('原始信号频谱');
 
grid on;axis([0 6000 0 400])
 
subplot(2,1,2)
 
plot(f2,abs(Y2));
 
xlabel('频率/Hz');ylabel('幅度');
 
title('过滤后信号的频谱');
 
grid on;axis([0 6000 0 100])
 


 
7.巴特沃斯低通滤波时频分析
 


 
clear
 
Fs=44100;
 
[s1,Fs,bits]=wavread('D:\222.wav');
 
x=s1(:,1);x1=x(40001:64576)
 
l=length(x1);
 
N=6;
 
m=downsample(x1,N);%降抽样后 Fs=7350hz
 
z=length(m);%4096
 
[tfr, t, f] = tfrstft(m);
 
figure(13)
 
contour(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
 
xlabel('时间t');
 
ylabel('频率f');
 
title('等高线图（滤波前）');
 
grid on;
 
figure(14)
 
mesh(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
 
xlabel('时间t');
 
ylabel('频率f');
 
zlabel('幅值A');
 
title('三维图（滤波前）');
 
grid on;
 
fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
 
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
 
[n,wn]=ellipord(wp,ws,rp,rs,'s');
 
[b,a]=ellip(n,rp,rs,wn,'s');
 
[num,den]=bilinear(b,a,Fs);
 
[h,w]=freqz(num,den,512,Fs);
 
y=filter(num,den,x);
 
x1=y(40001:64576)
 
l=length(x1);
 
N=6;
 
m=downsample(x1,N);%降抽样后 Fs=7350hz
 
z=length(m);%4096
 
[tfr,t,f] = tfrstft(m);
 
figure(15)
 
contour(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
 
xlabel('时间t');
 
ylabel('频率f');
 
title('等高线图（滤波后）');
 
grid on;
 
figure(16)
 
mesh(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
 
xlabel('时间t');
 
ylabel('频率f');
 
zlabel('幅值A');
 
title('三维图（滤波后）');
 
grid on;
 

 

 

Matlab数字滤波器的设计（IIR与FIR)（内含Matlab代码）
 
FIR滤波器的优缺点
 优点：相位线性度好，处理速度快，没有反馈回路稳定性强于IIR。
 缺点：FIR幅频特性精度较之于IIR低
 IIR滤波器
 优点：相同阶数下IIR滤波效果更好
 缺点：相位非线性程度较高，矫正时需要双向滤波进行矫正，不易控制。
 FIR滤波器最关键的在于窗函数：这里介绍
 汉明窗：
 FIR——Hamming窗
 
较低的旁瓣幅度，尤其是第一旁瓣；
旁瓣幅度要下降得快，以利于增加阻带衰减；
主瓣宽度要窄，这样滤波器过渡带较窄。
 N=40
 泄露系数0.04%，相对旁瓣抑制-42.1dB，主瓣宽度0.0625 （归一化后的频率）
 
 N越大，计算时间越长，群时延越长。取N=40（根据不同的情况，滤波器阶数可以取不同的值，这里我们对一段特定的音频信号进行滤波，因此取40）
 N越大，计算时间越长，群时延越长。取N=40
 

 滤波结果展示：
 
 蓝色为滤波前信号；红色为FIR滤波后的信号
 利用IIR椭圆滤波器对音频信号进行滤波：
 
 滤波结果展示：
 直接对歌曲进行IIR_ Eliptic（椭圆）滤波
 FIR滤波（黄色），IIR滤波后（紫色）
 对于不同的音频信号，对原始频谱进行分析后，只需要改动滤波器的设计参数，即可实现数字滤波。
 下面请看代码：
 
[x,fs]=audioread('C:\Users\a\Desktop\大三上\精密仪器\讨论课题目\新闻(加噪声).wav');
n=length(x);
  y=fft(x,n);       %做FFT变换
  f=fs*(0:n/2-1)/n;
  subplot(2,1,1);
  plot(x);
 axis([0 1400000 -1.5 1.5]);
  title('原始新闻信号采样后的时域波形'); 
xlabel('时间轴')
ylabel('幅值A')
subplot(2,1,2); 
plot(f,abs(y(1:n/2))); 
axis([0 1000 0 5000]);
title('原始新闻信号采样后的频谱图'); 
xlabel('频率Hz'); 
ylabel('频率幅值');
% n=0:N-1;
% w=2*n*pi/N;
FIRsignal=IIR_BANDPASSCheby;
outsignal=filter(FIRsignal,x);
%  IIRsignal=IIR_BANDPASSEliptic;
%  outsignal2=filter(IIRsignal,outsignal);
% set(gca,'fontsize',15,'tickdir','out')
% axis([0 5.0E-4 -2 2])
% xlabel('Time');ylabel('Amplitude')
% legend('FIR低通数字滤波器后的信号')  
 m=fft(outsignal,n);
 f=48000*(0:n/2-1)/n;
 figure
 subplot(211);
  plot(outsignal);
   axis([0 1400000 -1.5 1.5]);
  title('新闻滤波信号时域波形');
xlabel('时间轴')
ylabel('幅值A')
 subplot(212);
 plot(f,abs(m(1:n/2)))
 axis([200 2700 0 5000]);
  title('新闻滤波信号频谱图'); 
xlabel('频率Hz'); 
ylabel('频率幅值');
    sound(1.5*outsignal,44100);
%    audiowrite('C:\Users\a\Desktop\精密仪器\讨论课题目\新闻(IIR滤波后).wav',outsignal2,44100);
%    sound(x,44100);
 
第一次发，欢迎大家指正交流，不喜勿喷！！！！