1、1课 程 设 计课程名称_数字信号处理_题目名称_应用 Matlab 对语音信_号进行频谱分析及滤波学生学院_信息工程学院_专业班级_xx 级通信工程 xx 班_学 号_ 学生姓名_ 指导教师_ 2010 年 01 月 04 日2题目名称 应用 Matlab 对语音信号进行频谱分析及滤波学生学院 信息工程学院专业班级 xx 级通信 xx 班姓 名 xxxxxxxx学 号 xxxxxxxx一、课程设计目的数字信号处理是一门以算法为核心，理论和实践性较强的学科。是电子信息工程、通信工程专业、电子信息科学与技术专业的一门重要的专业技术基础课。数字信号处理课程是在学习完数字信号处理的相关理论后，进行的
2、综合性训练课程，其目的是：1、 使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；2、 增强学生应用 Matlab 语言编写数字信号处理的应用程序及分析、解决实际问题的能力；三、课程设计内容录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。下面对各步骤加以具体说明。21 语音信号的采集利用 Wi
3、ndows 下的录音机，录制一段自己的话音，时间在 1 s 内。然后在Matlab 软件平台下，利用函数 wavread 对语音信号进行采样，记住采样频率和采样点数。通过 wavread 函数的使用，学生很快理解了采样频率、采样位数等概念。22 语音信号的频谱分析首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。23 设计数字滤波器和画出其频率响应给出各滤波器的性能指标：3(1)低通滤波器性能指标 fb1 000 Hz，fc 1 200 Hz，As100 dB，Ap1 dB。(2)高通滤波器性能指标 fc4 800 Hz，fb5 000
4、Hz As100 dB，Ap1 dB。(3)带通滤波器性能指标 fb11 200 Hz，fb23 000 Hz，fc1 1 000 Hz， fc23 200 Hz，As100 dB，Ap1 dB。用窗函数法和双线性变换法设计上面要求的 3 种滤波器。在 Matlab 中，可以利用函数 fir1 设计 FIR 滤波器，可以利用函数 butte，cheby1 和 ellip 设计 IIR滤波器；利用 Matlab 中的函数 freqz 画出各滤波器的频率响应。24 用滤波器对信号进行滤波用自己设计的各滤波器分别对采集的信号进行滤波，在 Matlab 中，FIR 滤波器利用函数 fftfilt 对信
5、号进行滤波，IIR 滤波器利用函数 filter 对信号进行滤波。25 比较滤波前后语音信号的波形及频谱要求学生在一个窗口同时画出滤波前后的波形及频谱。26 回放语音信号在 Matlab 中，函数 sound 可以对声音进行回放。其调用格式：sound(x，fs，bits)；可以感觉滤波前后的声音有变化。27 设计系统界面设计处理系统的用户界面。在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。四、设计结果分析双线性变换法分析：双线性变换法的优点：不会有高于折叠频率的分量，避免了出现频率响应混叠响应。双线性变换法的缺点：1 一个线性相位的模拟滤波器经
6、双线性变换后就得到非线性相位的数字滤波器，不再保持原有的线性相位。2、这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，几某一频率段的幅频响应近似等于某一常数，否则变换所产生的数字滤波器幅频响应相当于原模拟滤波器的幅频响应会有畸变。窗函数法分析：窗函数法的优点：有闭合形式的公式可循，设计实用简单。窗函数法的缺点：通带、阻带的截止频率不易控制。4结果总结:只有选择的 fc、fb 与抽样频率 Fs 相匹配，才能得到与原声相近的声音输出。567891011121314五、设计心得体会在做这个课程设计之前，我就已接触过 MATLAB，学习过它的一些基本操作。但这次课程设计中涉及到很多 MATL
7、AB 的一些专门函数，对之前很少使用 MATLAB 的我来说，是个不小的挑战。在不久前的数字信号处理实验中，我又学习到一些用MATLAB 设计数字滤波器的方法，知道了一些专门函数。但即便如此，刚开始设计时还是觉得无处下手。为此我又重新温习了一遍教材中有关 IIR、FIR 滤波器的内容，自己动手计算了一次实验手册上的 IIR 低通滤波器的有关参数。在深感滤波器参数计算繁琐的情况下，我又上网搜索了一些用 MATLAB 设计数字滤波器的例程。只是例程毕竟是别人写的，在没多少注解的情况下，要读懂并不容易。为此，我再专门对一些感到疑惑的函数搜索一番，弄明白他们的使用形式。15最后，在众多例程中，我挑选了
8、一个接近该课程设计要求的 IIR 低通滤波器例程作参考。对其中的函数参数进行修改，得到了一个完全符合课程设计要求的 IIR 低通滤波器源程序。再依葫芦画瓢，接连写出其它的滤波器源程序。最终，我参考实验手册中例程框架，把已写出的几个滤波器源程序串在一起，得到了成品。经过这次课程设计，我对 IIR、FIR 数字滤波器的设计步骤加深了认识。同时，我也对 MATLAB 软件的强大感到吃惊。在数学研究方面，MATLAB 的确是首屈一指。有 MATLAB 的强大函数库支持，原本由人手计算的诸多滤波器参数可交由电脑处理，节省了大量时间。对 MATLAB 功能认识的进一步加深，我在今后的学习与工作中也能提高自
9、己的效率。六、参考资料1、 数字信号处理及 MATLAB 实现 余成波 杨如民等编著 清华大学出版社出版2、 数字信号处理教程 程佩青 清华大学出版社出版3.、 数字信号处理基础及 MATLAB 实现 周辉 董正宏 编著 北京希望电子出版社出版3、 MATLAB 在数字信号处理中的应用 薛年喜 主编 清华大学出版社七、源程序代码clear;close all;z1,fs,bits=wavread(shiyan.wav)%读取存放于同一文件夹中的声音文件length=input(length=);%D 是对原声取样长度y1=z1(1:length);%y1 是原声时域波形Y1=fft(y1);%
10、Y1 是原声 FFT 变换i=1;K=menu(请选择滤波器类型,IIR 低通滤波器,IIR 高通滤波器,.IIR 带通滤波器,FIR 低通滤波器,FIR 高通滤波器,.FIR 带通滤波器,原声,IIR 带阻滤波器, FIR 带阻滤波器,Exit);while(K=10)16if K=1%IIR 低通滤波器fb=input(fb=);fc=input(fc=);As=input(As=);Ap=input(Ap=);Fs=input(Fs=);wc=2*fc/Fs;wb=2*fb/Fs;N,ws=buttord(wb,wc,Ap,As)b,a=butter(N,ws)figure(i);i=i
11、+1;freqz(b,a,512,Fs);title(IIR 低通滤波器幅频相频图);x=filter(b,a,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(IIR 低通滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(IIR 低通滤波后信号波形);sound(x,fs,bits);elseif K=2%IIR 高通
12、滤波器fc=input(fc=);fb=input(fb=);As=input(As=);Ap=input(Ap=);Fs=input(Fs=);17wc=2*fc/Fs;wb=2*fb/Fs;N,ws=buttord(wc,wb,Ap,As);b,a=butter(N,ws,high);figure(i);i=i+1;freqz(b,a,512,Fs);title(IIR 高通滤波器幅频相频图);x=filter(b,a,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot
13、(2,2,2);plot(abs(X);title(IIR 高通滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(IIR 高通滤波后信号波形);sound(x,fs,bits);elseif K=3%IIR 带通滤波器fb1=input(fb1=);fb2=input(fb2=);fc1=input(fc1=);fc2=input(fc2=);As=input(As=);Ap=input(Ap=); Fs=input(Fs=);wc=2*fc1/Fs,2*fc2/Fs;wb=2*fb1/Fs
14、,2*fb2/Fs;N,ws=buttord(wb,wc,Ap,As);b,a=butter(N,ws,bandpass);figure(i);i=i+1;18freqz(b,a,512,Fs);title(IIR 带通滤波器幅频相频图);x=filter(b,a,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(IIR 带通滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);s
15、ubplot(2,2,4);plot(x);title(IIR 带通滤波后信号波形);sound(x,fs,bits);elseif K=4%FIR 低通滤波器fb=input(fb=);fc=input(fc=);As=input(As=);Ap=input(Ap=);Fs=input(Fs=);wc=2*pi*fc/Fs;wb=2*pi*fb/Fs;wdel=wc-wb;%确定过渡带宽 wdelbeta=0.1102*(As-8.7);%求凯泽窗参数 betaN=ceil(As-7.95)/2.286/wdel);%求凯泽窗阶数 Nwn=kaiser(N,beta);%凯泽窗函数ws=(w
16、b+wc)/2/pi;B=fir1(N-1,ws,wn);figure(i);i=i+1;freqz(B,1);title(FIR 低通滤波器幅频相频特性);x=fftfilt(b,z1);19X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(FIR 低通滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(FIR 低通滤波后信号波
17、形);sound(x,fs,bits);elseif K=5%FIR 高通滤波器fc=input(fc=);fb=input(fb=);As=input(As=);Ap=input(Ap=);Fs=input(Fs=);wc=2*pi*fc/Fs;wb=2*pi*fb/Fs;wdel=abs(wc-wb);%确定过渡带宽 wdelbeta=0.1102*(As-8.7);%求凯泽窗参数 betaN=ceil(As-7.95)/2.286/wdel);%求凯泽窗阶数 Nwn=kaiser(N,beta)%凯泽窗函数ws=(wb+wc)/2/pi;B=fir1(N-1,ws,high,wn);fi
18、gure(i);i=i+1;freqz(B,1);title(FIR 高通滤波器幅频相频特性);x=fftfilt(B,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);20title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(FIR 高通滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(FIR 高通滤波后信号波形);sound(x,fs,bits);elseif K=6%
19、FIR 带通滤波器fb1=input(fb1=);fb2=input(fb2=);fc1=input(fc1=);fc2=input(fc2=);As=input(As=);Ap=input(Ap=); Fs=input(Fs=);wb1=2*pi*fb1/Fs;wc1=2*pi*fc1/Fs;wb2=2*pi*fb2/Fs;wc2=2*pi*fc2/Fs;wdel=wb1-wc1;%确定过渡带宽 wdelbeta=0.1102*(As-8.7);%求凯泽窗参数 betaN=ceil(As-7.95)/2.286/wdel);%求凯泽窗阶数 Nwn= kaiser(N,beta);%凯泽窗函数
20、ws =(wb1+wc1)/2/pi,(wb2+wc2)/2/pi;B=fir1(N-1,ws,wn);figure(i);i=i+1;freqz(B,1);title(FIR 带通滤波器幅频相频特性);x=fftfilt(B,z1);X=fft(x,length);figure(i);i=i+1;21subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(FIR 带通滤波后信号频谱)subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot
21、(x);title(FIR 带通滤波后信号波形);sound(x,fs,bits);elseif K=7%原声figure(i);subplot(2,1,1);plot(abs(Y1);title(原始语音信号 FFT 频谱)subplot(2,1,2);plot(y1);title(原始语音信号时域波形)xlabel(时间 n);ylabel(幅值);sound(z1,fs,bits);elseif K=8%IIR 带阻滤波器fb1=input(fb1=);fb2=input(fb2=);fc1=input(fc1=);fc2=input(fc2=);As=input(As=);Ap=inp
22、ut(Ap=); Fs=input(Fs=);wc=2*fc1/Fs,2*fc2/Fs;wb=2*fb1/Fs,2*fb2/Fs;22N,ws=buttord(wb,wc,Ap,As);b,a=butter(N,ws,stop);figure(i);i=i+1;freqz(b,a,512,Fs);title(IIR 带阻滤波器幅频相频图);x=filter(b,a,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(IIR
23、带阻滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(IIR 带阻滤波后信号波形);sound(x,fs,bits);elseif K=9%FIR 带阻滤波器fb1=input(fb1=);fb2=input(fb2=);fc1=input(fc1=);fc2=input(fc2=);As=input(As=);Ap=input(Ap=); Fs=input(Fs=);wb1=2*pi*fb1/Fs;wc1=2*pi*fc1/Fs;wb2=2*pi*fb2/Fs;wc2=2*pi*fc2/F
24、s;wdel=wb1-wc1;%确定过渡带宽 wdelbeta=0.1102*(As-8.7);%求凯泽窗参数 betaN=ceil(As-7.95)/2.286/wdel);%求凯泽窗阶数 N23wn= kaiser(N,beta);%凯泽窗函数ws =(wb1+wc1)/2/pi,(wb2+wc2)/2/pi;B=fir1(N-1,ws,stop,wn);figure(i);i=i+1;freqz(B,1);title(FIR 带阻滤波器幅频相频特性);x=fftfilt(B,z1);X=fft(x,length);figure(i);i=i+1;subplot(2,2,1);plot(abs(Y1);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);title(FIR 带阻滤波后信号频谱)subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(FIR 带阻滤波后信号波形);sound(x,fs,bits); endK=menu(请选择滤波器类型,IIR 低通滤波器,IIR 高通滤波器,.IIR 带通滤波器,FIR 低通滤波器,FIR 高通滤波器,.FIR 带通滤波器,原声,IIR 带阻滤波器, FIR 带阻滤波器,Exit);endclose all;
