语音信号处理

1、语音信号的数字滤波处理 摘要本课程设计介绍了基于Matlab的对语音信号采集、处理及滤波器的设计，并使之实现的过程。运用课程中的基本概念、基本原理、基本分析方法，用Matlab进行数字语音信号处理，并阐述了课程设计的具体方法、步骤和内容。综合运用本课程的理论知识进行频谱分析以及滤波器设计，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。此外，系统实现了对数字信号处理的工作原理实施设计方法；用双线性变换法设计 IIR 数字滤波器，利用数字滤波器对信号进行滤波，对数字滤波器进行计算机仿真方法，并对设计结果加以分析。根据有关的频谱特性，直接法设计FIR数字滤波器，选择窗函

2、数，设计带通、高通、低通滤波方法进行滤波。 关键词： IIR滤波器；FIR滤波器；切比雪夫滤波器；汉宁窗；目录1绪论2设计原理2.1卷积运算2.2采样定理2.3用双线性变换法设计IIR数字滤波器2.4用窗函数法设计FIR滤波器3设计的前期准备.3.1语音信号的采集3.2线性卷积3.3采样定理4滤波器设计流程图5数字滤波器设计及其应用5.1含噪声语音信号的频谱分析5.2双线性变换法设计IIR滤波器5.3窗函数设计FIR滤波器6心得与体会7参考文献1绪论      数字滤波器可以在语音信号分析中对声音进行处理，可以滤出不要的噪声

3、，使声音更加清楚。本设计通过对语音信号进行采集，对语音信号进行时域与频域的分析，然后给语音信号加上噪声，通过切比雪夫滤波器进行高通、低通、带通的滤波。通过汉宁窗对声音进行过滤。然后对声音进行回放，对比前后声音信号的差异。实现滤波功能。理论依据：根据设计要求分析系统功能，掌握设计中所需理论（采样频率、采样位数的概念，采样定理； 时域信号的FFT分析；数字滤波器设计原理和方法，各种不同类型滤波器的性能比较），阐明设计原理。信号采集：采集语音信号，并对其进行FFT频谱分析，画出信号的时域波形图和频谱图。构造受干扰信号并对其进行FFT频谱分析：对所采集的语音信号加入干扰噪声，对语音信号进行回

4、放，感觉加噪前后声音的变化，分析原因，得出结论。并对其进行FFT频谱分析，比较加噪前后语音信号的波形及频谱，对所得结果进行分析，阐明原因，得出结论。 数字滤波器设计：根据待处理信号特点，设计合适数字滤波器，绘制所设计滤波器的幅频和相频特性。信号处理：用所设计的滤波器对含噪语音信号进行滤波。对滤波后的语音信号进行FFT频谱分析。画出处理过程中所得各种波形及频谱图。对语音信号进行回放，感觉滤波前后声音的变化。比较滤波前后语音信号的波形及频谱，对所得结果和滤波器性能进行频谱分析，阐明原因，得出结论。2设计原理2.1卷积运算 卷积和乘积运算在频域和时域是一一对应的，两个信号在时域的卷积可以转

5、化为求两者在频域的乘积后再反变换，同理在频域的卷积等时域的乘积。而信号的频域求解有快速傅里叶FFT算法。卷积与傅里叶变换有着密切的关系。利用一点性质，即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，能使傅里叶分析中许多问题的处理得到简化。由卷积得到的函数f*g 一般要比f 和g 都光滑。特别当g 为具有紧支集的光滑函数，f 为局部可积时，它们的卷积f * g 也是光滑函数。利用这一性质，对于任意的可积函数f，都可以简单地构造出一列逼近于f 的光滑函数列，这种方法称为函数的光滑化或正则化。卷积的概念还可以推广到数列、测度以及广义函数上去。2.2采样定理采样定理，又称香农采样定理，奈奎斯特采

6、样定理，是信息论，特别是通讯与信号处理学科中的一个重要基本结论.E.T.Whittaker(1915年发表的统计理论),克劳德·香农与Harry Nyquist都对它作出了重要贡献。另外,V.A.Kotelnikov也对这个定理做了重要贡献。采样是将一个信号（即时间或空间上的连续函数）转换成一个数值序列（即时间或空间上的离散函数）。采样定理指出，如果信号是带限的，并且采样频率高于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。带限信号变换的快慢受到它的最高频率分量的限制，也就是说它的离散时刻采样表现信号细节的能力是有限的。采样定理是指，如果信号带宽不到采样频率的一半

7、（即奈奎斯特频率），那么此时这些离散的采样点能够完全表示原信号。高于或处于奈奎斯特频率的频率分量会导致混叠现象。大多数应用都要求避免混叠，混叠问题的严重程度与这些混叠频率分量的相对强度有关。 2.3用双线性变换法设计IIR数字滤波器 脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=esT1将此横带变换到整个

8、Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。S平面上与Z平面的成非线性的正切关系。在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。但是双线性变换的这个特点是靠频率的严重非线性关系

9、而得到的，由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变。对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正

10、好映射到所需要的数字频率上。 2.4用窗函数法设计FIR滤波器 根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为w，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数

11、wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样N点，采用离散傅里叶反变换（IDFT）即可求出。根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程

12、序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。3设计前期准备 3.1语音信号的采集1.语音信号的采集：利用windows下的c盘（C:WINDOWSMedia recycle，16位，双声道），然后将音频文件保存“H:recycle .wav”2在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。3语音信号的频谱分析首先画出语音信号的时域波形z1=wavread(' H:recycle .wav ');plot(z1);title('原始信号时域波形');xlabel('数据

13、采集单位 n');ylabel('amplitude');图像输出如图1.1图1.1语音信号时域波形语音信号的频谱特性图1.2z1=wavread('H:recycle.wav');plot(z1);xlabel('Frequency（hz）');ylabel('Amplitude');title('语音信号频谱特性'); 图1.2语音信号频谱特性 3.2线性卷积   任意两个序列x1(n)、x2(n)，指定x1(n)为自己的学号,例如x1(n)=2,0,0

14、,7,8,4,2,5,0,1,2,3。x2(n)的内容和长度自选进行线性卷积。   线性系统的含义，就是，这个所谓的系统，带来的输出信号与输入信号的数学关系式之间是线性的运算关系。因此，实际上，都是要根据我们需要待处理的信号形式，来设计所谓的系统传递函数，那么这个系统的传递函数和输入信号，在数学上的形式就是所谓的卷积关系。卷积关系最重要的一种情况，就是在信号与线性系统或数字信号处理中的卷积定理。利用该定理，可以将时间域或空间域中的卷积运算等价为频率域的相乘运算，从而利用FFT等快速算法，实现有效的计算，节省运算代价。线性卷积程序如下clearxls=2,0,0,7

15、,8,4,2,5,0,1,4,1;%输入X数组lx=length(xls);hls=1,1,1,1;%输入h数组lh=length(hls);lmax=max(lx,lh);if lx>lh    nx=0;nh=lx-lh;elseif  lx<lh  nx=lh-lx;nh=0;else  nx=0;nh=0;enddt=0.5;%运算时间间隔lt=lmax;%取长者为补0长度为基准x=zeros(1,lt),xls,zeros(1,nx),zeros(1,lt);%将X补锝与h

16、同长，在两边补上同长度的0tl=(-lt+1:2*lt)*dt;h=zeros(1,2*lt),hls,zeros(1,nh);hf=fliplr(h);%将h的左右翻转，记为hfy=zeros(1,3*lt);%y数组初始化for k=0:2*lt;    p=zeros(1,k),hf(1:end-k);%hf向又循环移位    y1=x.*p*dt;%输入和翻转移位脉冲函数相乘，在乘dt    yk=sum(y1);%相加，等于积分   &

17、#160;y(k+lt+1)=yk;%将结果放入数组Y    subplot(4,1,1);stem(tl,x)%用stem表示每次卷积结果       axis(-lt*dt,2*lt*dt,min(x),max(x),hold on;ylabel('x(t)');     title('卷积');    subplot(4,1,2);stem(tl,p)

18、0;   axis(-lt*dt,2*lt*dt,min(p),max(p),    ylabel('h(k-t)');    subplot(4,1,3);stem(tl,y1)    axis(-lt*dt,2*lt*dt,min(y1),max(y1)+eps),    ylabel('s=x*h(k-t)');    subpl

19、ot(4,1,4);stem(k*dt,yk)    axis(-lt*dt,2*lt*dt,floor(min(y)+eps),ceil(max(y)+eps),hold on;    ylabel('y(k)=sum(s)*dt)');    if k=round(0.8*lt)   disp('pause press any key to countinue'),pause   

20、 else pause(1),    endend线性卷积输出如图3.2图3.2线性卷积 3.3采样定理采样位数即采样值或取样值，用来衡量声音波动变化的参数，采样频率是指在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精。采样程序如下： % 时域采样理论验证Tp=64/1000;

21、%产生M长采样序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=41;alph=41*20.5*pi;omega=pi*41*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);  %M点FFTxnt)yn='xa(nT)'subplot(3,2,1);tstem(xnt,yn);box on;title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,a

22、bs(Xk);title('(a) T*FTxa(nT),Fs=1000Hz');xlabel('f(Hz)');ylabel('幅度');axis(0,180,0,1.2*max(abs(Xk)% Fs=300;T=1/Fs;Fs=300;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);  %M点FFTxnt)yn='xa(nT)'subplot(3,2,3);tstem(xnt,yn);box

23、 on;title('(b) Fs=300Hz');k=0:M-1;fk=k/Tp;subplot(3,2,4);plot(fk,abs(Xk);title('(b) T*FTxa(nT),Fs=300Hz');xlabel('f(Hz)');ylabel('幅度');axis(0,180,0,1.2*max(abs(Xk)% Fs=200;T=1/Fs;Fs=200;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);

24、  %M点FFTxnt)yn='xa(nT)'subplot(3,2,5);tstem(xnt,yn);box on;title('(c) Fs=200Hz');k=0:M-1;fk=k/Tp;subplot(3,2,6);plot(fk,abs(Xk);title('(c) T*FTxa(nT),Fs=200Hz');xlabel('f(Hz)');ylabel('幅度');axis(0,180,0,1.2*max(abs(Xk)时域采样理论的验证程序运行如图3.3所示。由图可见，采样序列的频

25、谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。   图3.3时域采样4滤波器设计流程图 本课程设计主要是对语音信号进行加噪声，滤波。首先要对声音进行采集，才可以进行信号处理。有了信号接下来要考虑的就是声音进行加噪声。可以插入一个白噪声或者脉冲信号。观察声音在频谱和时域在加噪声前后的变化，及加噪声与没加噪声的区别。可以通过对语音信号的时域波形，频谱进行比较。然后对语音信号进行频谱分

26、析，对比。再通过MATLAB，根据有关的频谱特性，采用间接法设计IIR数字滤波器切比雪夫滤波器进行低通、带通、高通滤波，看噪声过滤情况。还可以根据有关的频谱特性，采用直接法设计FIR数字滤波器汉宁窗进行低通、带通、高通滤波，对比滤波情况。然后对滤波后波形进行频谱分析、最后在Matlab中,函数sound可以对声音进行回放。其调用格式：sound(x,fs,bits);可以感觉滤波前后的声音有变化。录制或打开一个WAV音频文件开始窗函数设计FIR滤波器用双线性变换法设计IIR滤波器频谱分析录音回放对比结束加入噪声信号对信号进行频谱时域分析低涌滤波高通滤波带通滤波  5数字滤波

27、器设计及其应用 5.1污染信号的频谱分析对信号进行采集后，往往有噪声。对于这些噪声应该怎么处理。我们可以用滤波器对信号进行滤波、对频谱进行分析。通过matlab的绘图功能进行比较，观察加噪声与不加噪声频谱的不同。1、污染信号的频谱分析程序如下Clearz1,fs,bits=wavread('H:recycle');y1=z1(1:8192);t=0:1/22050:8191*(1/22050);z=sin(20000*pi*t);m=y1+0.00000005*z'u=fft(m,8192);f=22050*(0:4096)/8192;wp=5000/1102

28、5;ws=6000/11025;n,wn=buttord(wp,ws,3,20);b,a=butter(n,wn);freqz(b,a,8192,22050);s=filter(b,a,m);subplot(2,2,2);plot(f,abs(u(1:4097);title('滤波前信号频谱 ');xlabel('Frequency（hz）');ylabel('Amplitude')subplot(2,2,1);plot(m(1:8192);title('滤波前信号 '); xlabel('n');ylabel(&

29、#39;Amplitude') subplot(2,2,3);plot(s(1:8192);title('滤波后信号');y=fft(s,8192);%对x进行256点的fftsubplot(2,2,4);plot(f,abs(y(1:4097);%画出频域内的信号title('滤波后信号频谱 '); xlabel('Frequency（hz）');ylabel('Amplitude') 2、加噪声语音信号时域波形和频谱如图4.1 图5.1加噪声语音信号时域波形与频谱 5.

30、2双线性变换法设计低通滤波器      Chebyshev 滤波器则比Butterworth 滤波器的截至特性要好，在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。但通带处的幅值有振荡。对于数字滤波器而言，可以采用不同阶数逼近相应滤波器，滤波器性能还与滤波器的阶数有关，一般而言，阶数越高，则逼近越精确，但计算代价也随之上升，所以性能与代价总需要寻求一个平衡点。本设计用Chebyshev 滤波器。4.2.1选用cheby1高通在设计数字滤波器时，cheb

31、y1函数中的参数Wn是一个相对量，其定义区间为Wn0,1，其中1对应于0.5fs，fs为采样频率（单位Hz）；在设计滤波器时，Wn采用真实频率，单位为Hz。 b,a = cheby1(n,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶b,a = cheby1(n,Rp,Wn,'ftype')用于设计数字高通和数字带阻滤波器，即ftype = high时，返回截止频率为Wn的数字高通滤波器； 程序设计如下： %高通滤波器选用cheby1clear;close allz1,fs,bits=wavread('H:recycle.wav')y

32、1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=22050;wc=2*fc/Fs;wb=2*fp/Fs;% 设置指标参数n,wp=cheb1ord(wc,wb,Ap,As);%计算切比雪夫1型模拟高通滤波器阶数和同%带边界频率b,a=cheby1(n,Ap,wp,'high');% 切比雪夫1型滤波器系统函数系数%作图部分figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0

33、,4000,0,2.4);title('滤波前信号频谱');xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,2);plot(abs(X);axis(0,4000,0,1.2);title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,3);plot(z1);title('滤波前信号波形'); xlabel('n&#

34、39;);ylabel('Amplitude')subplot(2,2,4);plot(x);title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')sound(x,fs,bits); 图形分析如下图5.2.1.1 图5.2.1.2                     

35、   图5.2.1.1                       图5.2.1.2 cheby1高通滤波 5.2.2选用cheby1低通cheby1函数中的参数Wn是一个相对量，其定义区间为Wn0,1，其中1对应于0.5fs，fs为采样频率（单位Hz）；在设计滤波器时，Wn采用真实频率，单位为Hz。b,a = cheby1(n

36、,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶ChebyshevI型数字低通滤波器，通带内波纹为Rp。b、a分别为滤波器传递函数的分子和分母系数向量（降幂排列）。程序设计如下：%双线性变换法设计低通滤波器%选用cheby1clear;close allz1,fs,bits=wavread('H:recycle.wav')y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=40;Ap=1; ;Fs=22050; % 设置指标参数wc=2*fc/Fs;wb=2*fp/Fs; % 设置指标参数n,wp=cheb1ord(wc,wb,Ap,

37、As); %计算切比雪夫1型模拟低通滤波器阶数和同%带边界频率 b,a=cheby1(n,Ap,wp，low); % 切比雪夫1型滤波器系统函数系数figure(1);h,w=freqz(b,a);plot(w,abs(h);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,2.0);title('滤波前信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');sub

38、plot(2,2,2);plot(abs(X);axis(0,1000,0,2.0);title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,3);plot(z1);title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')subplot(2,2,4);plot(x);title('滤波后号波形'); xlabel('n');ylabel

39、('Amplitude')sound(x,fs,bits); 图形分析如下图5.2.2.1 图5.2.2.2                            图5.2.2.1图5.2.1.2 cheby1低通滤波 5.2.2选用cheby1带通cheby1函数中的参数Wn采用真

40、实频率，单位为Hz。b,a = cheby1(n,Rp,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶ChebyshevI型数字低通滤波器，当Wn为二元向量，即Wn=W1 W2 (W1<W2)时，b,a = cheby1(n,Rp,Wn)返回一个2n阶数字带通滤波器，其通带为W1<< W2。 程序设计如下：%带通滤波器选用cheby1clear;close allz1,fs,bits=wavread('H:recycle.wav')y1=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=

41、3200 ;As=30;Ap=1; Fs=22050; wc=2*fc1/Fs,2* fc2/Fs;wb=2*fp1/Fs,2*fp2/Fs; % 设置指标参数 n,wpo=cheb1ord(wc,wb,Ap,As); %计算切比雪夫1型模拟带通滤波器阶数和同带边界频率 b,a=cheby1(n,Ap,wpo); % 切比雪夫1型滤波器系统函数系数%作图部分figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,3000,0,2

42、.5);title('滤波前信号频谱');xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,2);plot(abs(X);axis(0,3000,0,2.5);title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,3);plot(z1);title('滤波前信号波形'); xlabel('n');ylabel(&#

43、39;Amplitude')subplot(2,2,4);plot(x);title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')sound(x,fs,bits); 图形分析如下图5.2.2.1 图5.2.2.2图5.2.3.1图5.2.3.2 cheby1带通滤波 5.3窗函数设计滤波器汉宁窗又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是 3个 sine（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 /T，从而使旁瓣互相抵消，消去高频干扰和

44、漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。这里用的就是汉宁窗。5.3.1窗函数设计低通滤波器选择一个适当的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。根据窗函数最小阻带衰减的特性，可采用汉宁窗提供大于30dB的衰减。本设计选择汉宁窗，其过渡带为8/N，因此具有较小的阶次。       程序设计如下：clear;close allz1,fs,bits=wavread('H:recycle.wav')

45、y1=z1(1:8192);Y1=fft(y1);fp=1000 ;fc=1200  As=100 ;Ap=1 ;Fs=22050 ;wp=2*fp/Fs;ws=2*fc/Fs;DB=ws-wp;           %计算过渡带宽度N0=ceil(6.2*pi/DB); %根据表7.2.2汉宁窗计算所需h(n)长度N0N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数wc=(wp+ws)/2/pi; 

46、0;  %计算理想高通滤波器通带截止频率(关于归一化）hn=fir1(N-1,wc,'low',hanning(N);  %调用fir1计算高通FIRDFh(n)%以下是绘图部分figure(1);freqz(hn,1);x=fftfilt(hn,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,300,0,2.0);title('滤波前信号频谱');xlabel('Frequency（hz）');ylabel('Ampli

47、tude');subplot(2,2,2);plot(abs(X);axis(0,300,0,0.5);title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,3);plot(z1);title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')subplot(2,2,4);plot(x);title('滤波后信号波形'); xlabel('

48、n');ylabel('Amplitude') 图形分析如下图5.3.1.1 图5.3.1.2 图5.2.3.1                              图5.2.3.2 汉宁窗低通滤波器 5.3.2窗函数设计高通滤波器数字高通

49、滤波器的设计分两种情况，即无相移和相移为±/2两种情况。当无相移时，N为奇数时，所设计的FIR数字高通滤波器为型滤波器；当N为偶数时，为型滤波器。当相移为±/2时，N为奇数时，所设计的FIR数字高通滤波器为型滤波器。根据窗函数最小阻带衰减的特性，选择汉宁窗可达到 44dB的最小阻带衰减，其过渡带为6.2/N，因此具有较小的阶次。 程序设计如下：%窗函数设计高通滤波器clear;close allz1,fs,bits=wavread('H:recycle.wav')y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;

50、As=100;Ap=1; Fs=22050;wp=2*fp/Fs;ws=2*fc/Fs;DB=wp-ws;           %计算过渡带宽度N0=ceil(6.2*pi/DB); %根据表7.2.2汉宁窗计算所需h(n)长度N0N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数wc=(wp+ws)/2/pi;    %计算理想高通滤波器通带截止频率(关于归一化）hn=fir1(N-1,wc,

51、9;high',hanning(N);  %调用fir1计算高通FIRDFh(n)%以下是绘图部分figure(1);freqz(hn,1);x=fftfilt(hn,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,800,0,1.0);title('滤波前信号频谱');xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,2);plot(abs(X);axis(0,800,0,1.0

52、);title('滤波后信号频谱'); xlabel('Frequency（hz）');ylabel('Amplitude');subplot(2,2,3);plot(z1);title('滤波前信号波形'); xlabel('n');ylabel('Amplitude')subplot(2,2,4);plot(x);title('滤波后信号波形'); xlabel('n');ylabel('Amplitude')sound(x,fs,bits);&#

53、160;图形分析如下图5 图5.3.2汉宁窗低通滤波器 5.3.3窗函数设计带通滤波器 数字带通滤波器的设计也分两种情况，即无相移和相移为±/2两种情况。当无相移时，N为奇数时，所设计的FIR数字带通滤波器为型滤波器；当N为偶数时，为型滤波器。当相移为±/2时，N为奇数时，所设计的FIR数字带通滤波器为型滤波器；当N为偶数时，为型滤波器。汉宁窗通带纹波和阻带纹波均满足设计要求。这里选汉宁窗。程序设计如下：%窗函数设计带通滤波器clear;close allz1,fs,bits=wavread('H:recycle.wav')y1

54、=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000 ;fc1=1000  fc2=3200 ;As=100 ;Ap=1 ;Fs=22050 ;wp1=2*pi*fp1/Fs; wc1=2*pi*fc1/Fs; wp2=2*pi*fp2/Fs; wc2=2*pi*fc2/Fs;wdel=wp1-wc1;       %计算过渡带宽度N0=ceil(6.2*pi/wdel); %根据表7.2.2汉宁窗计算所需h(n)长度N0N=N0+mod(N0+1,2);   %确保h(n)长度N是奇数wc=(wp1+