应用Matlab对语音信号进行频谱分析及滤波
一、课程设计的内容
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。
二、课程设计的要求与数据
1、 学生能够根据设计内容积极主动查找相关资料; 2、 滤波器的性能指标可以根据实际情况作调整; 3、 对设计结果进行独立思考和分析; 4、 设计完成后,要提交相关的文档;
1) 课程设计报告书(纸质和电子版各一份,具体格式参照学校课程设计管理规定),
报告内容要涵盖设计过程、频谱图的分析.
2) 可运行的源程序代码(电子版)
5、 在基本要求的基础上,学生可以根据个人对该课程设计的理解,添加一些新的内
容;
6、 详细设计要求参照>课程设计指导手册.
三、课程设计应完成的工作
1、 语音信号的采集; 2、 语音信号的频谱分析; 3、 数字滤波器的设计; 4、 对语音信号进行滤波处理;
5、 对滤波前后的语音信号频谱进行对比分析;
四、课程设计进程安排
五、应收集的资料及主要参考文献
1、程佩青.数字信号处理教程.北京：清华大学出版社出版，2001
2、 刘方铭，姚震，韩国军等.数字信号处理实验指导书.广州：广东工业大学信息工程学院，2006
3、辅助信号处理技术与应用 编著：飞思科技产品研发中心 北京：电子工业出版社2005.3
发出任务书日期： 2007 年 12 月 21 日 指导教师签名：李学易
计划完成日期： 2007 年 12 月 28 日 基层教学单位责任人签章：
主管院长签章：
一、设计题目
应用Matlab对语音信号进行频谱分析及滤波。 二、设计目的
为了巩固所学的数字信号处理理论知识，使学生对信号的采集、处理、传输、显示和存储等有一个系统的掌握和理解，再者，加强学生对Matlab软件在信号分析和处理的运用。 三、设计内容
1 语音信号的采集 利用windows自带的录音机，录制一段语音，时间在1s以内。在matlab平台下，利用函数wavread对语音信号进行采样。通过wavread函数的使用，进一步理解采样频率，采样位数。
2 语音信号的频谱分析
首先画出由wavread函数采样后的信号的时域波形和它的频域响应波形；其次对语音信号进行FFT变换，得到FFT频谱特性曲线，与原语音信号的频谱特性曲线进行比较。 3 设计数字滤波器和画出其频率响应
用窗函数法和双线性变换法设计以下三种数字滤波器：
（1）低通滤波器性能指标fb=1000Hz，fc＝1200Hz， As＝100db，Ap＝1db。 （2）高通滤波器性能指标 fs＝4800Hz，fb＝5000Hz，As＝100db，Ap＝1db。
（3）带通滤波器性能指标 fb1＝1200Hz，fb2＝3000Hz，fc1＝1000Hz，fc2＝3200Hz，As＝100db，Ap＝1db。 4 系统界面设计
本系统界面主要采用if…else语句和menu菜单函数进行设计，如下图：
四、设计结果分析
1原始信号采样后时域图和频谱图： （1）原始语音信号时域图和频率响应图：
（2）原始语音信号频谱图及做FFT变换后频谱图：
2 IIR滤波器及通过IIR滤波器后信号的时域和频谱的变化情况并进行比较： （1）IIR低通滤波器：
分析：题目要求,最大衰减AP=1dB=0.89时fb=1000Hz，由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc=1200Hz，基本符合
分析：由于所取的采样点数比较大，滤波前后的频谱比较相近，但仔细对比仍然可以看出下图结果滤掉了高频成分。
（2）IIR高通滤波器：
分析：题目要求,最大衰减AP=1dB=0.89时fb=5000Hz，由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc=4800Hz，基本符合
分析：由图明显看出，该滤波器将低频成分滤掉了 （3）IIR带通滤波器：
分析：题目要求,最大衰减AP=1dB=0.89时fb1=1200Hz，fb2=3000Hz由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc1=1000Hz，fc2=3200Hz基本符合。
分析：上图很明显看出该滤波器将低频和高频成分滤掉了，只剩下通带里的频率。 3 FIR滤波器及通过FIR滤波器后信号的时域和频谱的变化并进行比较： （1）FIR低通滤波器：
（2）FIR高通滤波器：
（3）FIR带通滤波器：
五、设计心得体会
通过做应用Matlab对语音信号进行谱分析及滤波这个课程设计，加深了我对数字信号处理理论知识，尤其是采样频率，频谱特性和数字滤波器的设计；另一方面，使我在原来熟悉控制系统在Matlab运用的基础上进一步学习Matlab软件在信号分析和处理的运用。在整个课程设计过程中，从分析设计题目到设计程序框图再到具体的程序代码的编写，使我对课程设计有了比较熟悉的了解（这是我第一次做课程设计）。
分析课题题目：应用Matlab对语音信号进行谱分析及滤波，考虑到几个方面的内容：（1）加强对Matlab基础知识的学习，尤其是Matlab在信号处理方面的学习，这就需要多方面查找资料，如上互联网，图书馆等；（2）加强学习信号处理知识理论，尤其是FFT频谱变换理论和各种数字滤波器设计理论；（3）要突出重点，其重点是谱分析和数字滤波器的设计和分析结果。
明确设计步骤，设计程序框图，按照框图编写程序代码，调试修改及完善。在具体编写代码的过程中，遇到了不少问题，比如设计数字滤波器所需函数的使用，归一化与非归一化的问题，各种参数的确定等等，但是通过对资料仔细的通阅和不断的调试修改，终于将课程设计完成，自己个人觉得质量还算优秀：另外，我觉得我在设计的过程中还有一点是很重要的：就是要有对自己负责的心态，永不言弃。
六、参考文献
1、程佩青.数字信号处理教程.北京：清华大学出版社出版，2001
2、 刘方铭，姚震，韩国军等.数字信号处理实验指导书.广州：广东工业大学信息工程学院，2006
3、辅助信号处理技术与应用 编著：飞思科技产品研发中心 北京：电子工业出版社2005.3
七附录：
b=menu('请选择选项','原始信号采样后时域图和频谱图','FIR滤波器','IIR滤波器','退出'); if b==4 b==0; end if b==1
temp=menu('请选择选项','播放原始语音','原始语音时域图和频率响应图','原始语音频谱图及做FFT变换后频谱图','返回'); if temp==1 voice; main;
elseif temp==2 time; main; elseif temp==3 pingpu; main; else
main; end
elseif b==2
temp=menu('请选择选项','FIR低通滤波器','FIR高通滤波器','FIR带通滤波器','返回'); if temp==1 FIR_LP; main; elseif temp==2 FIR_HP; main; elseif temp==3 FIR_BP; main; else
main; end
elseif b==3
temp=menu('请选择选项','IIR低通滤波器','IIR高通滤波器','IIR带通滤波器','返回'); if temp==1
IIR_LP; main; elseif temp==2 IIR_HP; main; elseif temp==3 IIR_BP; main; else main end end
fs=22050;
x1=wavread('录音.wav');
wp1=2*pi*1200/fs;wp2=2*pi*3000/fs; ws1=2*pi*1000/fs;ws2=2*pi*3200/fs; Rp=1; Rs=100;
wp=(wp1+ws1)/2;ws=(wp2+ws2)/2; wdelta=wp1-ws1;
N=ceil(8*pi/wdelta); %取整 wn=[wp ws];
[b,a]=fir1(N,wn/pi,'bandpass'); figure(1)
freqz(b,a,512);
title('FIR带通滤波器'); f2=filter(b,a,x1); figure(2)
subplot(2,1,1) plot(x1)
title('FIR带通滤波器滤波前的时域波形'); subplot(2,1,2) plot(f2);
title('FIR带通滤波器滤波后的时域波形');
sound(f2,44100); %播放滤波后的语音信号 F0=fft(f2,1024); f=fs*(0:511)/1024; figure(3)
y2=fft(x1,1024); subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR带通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值'); subplot(2,1,2)
plot(f,abs(F0(1:512)));
title('FIR带通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');
fs=22050;
x1=wavread('录音.wav'); wp=2*pi*5000/fs; ws=2*pi*4800/fs; Rp=1; Rs=100;
wdelta=wp-ws;
N=ceil(8*pi/wdelta); %取整 wn=(wp+ws)/2;
[b,a]=fir1(N,wn/pi,'high'); figure(1)
freqz(b,a,512);
title('FIR高通滤波器'); f2=filter(b,a,x1); figure(2)
subplot(2,1,1) plot(x1)
title('FIR高通滤波器滤波前的时域波形'); subplot(2,1,2) plot(f2);
title('FIR高通滤波器滤波后的时域波形');
sound(f2,44100); %播放滤波后的语音信号 F0=fft(f2,1024); f=fs*(0:511)/1024; figure(3)
y2=fft(x1,1024); subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR高通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值'); subplot(2,1,2)
F2=plot(f,abs(F0(1:512)));
title('FIR高通滤波器滤波后的频谱')
xlabel('频率/Hz'); ylabel('幅值');
fs=10000;
x1=wavread('录音.wav'); wp=2*pi*1000/fs; ws=2*pi*1200/fs; Rp=1; Rs=100;
wdelta=ws-wp;
N=ceil(8*pi/wdelta); %取整 wn=(wp+ws)/2;
[b,a]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化截止频率 figure(1)
freqz(b,a,512);
title('FIR低通滤波器'); f2=filter(b,a,x1); figure(2)
subplot(2,1,1) plot(x1)
title('FIR低通滤波器滤波前的时域波形'); subplot(2,1,2) plot(f2);
title('FIR低通滤波器滤波后的时域波形');
sound(f2,44100); %播放滤波后的语音信号 F0=fft(f2,1024); f=fs*(0:511)/1024; figure(3)
y2=fft(x1,1024); subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR低通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值'); subplot(2,1,2)
F2=plot(f,abs(F0(1:512)));
title('FIR低通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');
fs=22050;
x1=wavread('录音.wav');
Fs=22050;
Ts=1/Fs;R1=30;
fb1=1200;fb2=3000;fc1=1000;fc2=3200;fs=22050; W1=2*fb1*pi/fs;W2=2*fc1*pi/fs; W3=2*fb2*pi/fs;W4=2*fc2*pi/fs; Wp=[W1,W3]; Ws=[W2,W4]; Rp=1; Rl=100;
Wp1=2/Ts*tan(Wp/2); %将模拟指标转换成数字指标 Ws1=2/Ts*tan(Ws/2);
[N,Wn]=cheb2ord(Wp1,Ws1,Rp,R1,'s'); %选择滤波器的最小阶数
[Z,P,K]=cheb2ap(N,Rl); %创建切比雪夫模拟滤波器 [Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2bp(Bap,Aap,2100*2*pi,1800*2*pi);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,W]=freqz(bz,az); %绘制频率响应曲线 figure(1)
plot(W*fs/(2*pi),abs(H)) grid
xlabel('频率／Hz') ylabel('频率响应幅度') title('IIR带通滤波器') f1=filter(bz,az,x1); figure(2)
subplot(2,1,1)
plot(x1) %画出滤波前的时域图 title('IIR带通滤波器滤波前的时域波形'); subplot(2,1,2)
plot(f1); %画出滤波后的时域图 title('IIR带通滤波器滤波后的时域波形');
sound(f1,44100); %播放滤波后的信号 F0=fft(f1,1024); f=fs*(0:511)/1024; figure(3)
y2=fft(x1,1024); subplot(2,1,1);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图 title('IIR带通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值'); subplot(2,1,2)
plot(f,abs(F0(1:512))); %画出滤波后的频谱图 title('IIR带通滤波器滤波后的频谱')
xlabel('频率/Hz');
ylabel('幅值');
fs=22050;
x1=wavread('录音.wav');
Fs=22050;
Ts=1/Fs;R1=50;
Wp=2*pi*5000/fs;
Ws=2*pi*4800/fs;
Rp=1;
Rl=100;
Wp1=2/Ts*tan(Wp/2); %将模拟指标转换成数字指标
Ws1=2/Ts*tan(Ws/2);
[N,Wn]=cheb2ord(Wp1,Ws1,Rp,Rl,'s'); %选择滤波器的最小阶数
[Z,P,K]=cheb2ap(N,Rl); %创建切比雪夫模拟滤波器
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2hp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz(bz,az); %绘制频率响应曲线
figure(1)
plot(W*fs/(2*pi),abs(H))
grid
xlabel('频率／Hz')
ylabel('频率响应幅度')
title('IIR高通滤波器')
f1=filter(bz,az,x1);
figure(2)
subplot(2,1,1)
plot(x1) %画出滤波前的时域图
title('IIR高通滤波器滤波前的时域波形');
subplot(2,1,2)
plot(f1); %画出滤波后的时域图
title('IIR高通滤波器滤波后的时域波形');
sound(f1,44100); %播放滤波后的信号
F0=fft(f1,1024);
f=fs*(0:511)/1024;
figure(3)
y2=fft(x1,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图
title('IIR高通滤波器滤波前的频谱')
xlabel('频率/Hz');
ylabel('幅值');
plot(f,abs(F0(1:512))); %画出滤波后的频谱图
title('IIR高通滤波器滤波后的频谱')
xlabel('频率/Hz');
ylabel('幅值');
fs=22050;
[x2]=wavread('录音.wav');
Ts=1/fs;R1=10;
wp=2*pi*1000/fs;
ws=2*pi*1200/fs;
Rp=1;
Rl=100;
wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标
ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,R1,'s'); %选择滤波器的最小阶数
[Z,P,K]=buttap(N); %创建butterworth模拟滤波器
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz(bz,az); %绘制频率响应曲线
figure(1)
plot(W*fs/(2*pi),abs(H))
grid
xlabel('频率／Hz')
ylabel('频率响应幅度')
title('IIR低通滤波器')
f1=filter(bz,az,x2);
figure(2)
subplot(2,1,1)
plot(x2) %画出滤波前的时域图
title('IIR低通滤波器滤波前的时域波形');
subplot(2,1,2)
plot(f1); %画出滤波后的时域图
title('IIR低通滤波器滤波后的时域波形');
sound(f1,44100); %播放滤波后的信号
F0=fft(f1,1024);
f=fs*(0:511)/1024;
figure(3)
y2=fft(x2,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图
title('IIR低通滤波器滤波前的频谱')
xlabel('频率/Hz');
ylabel('幅值');
F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图
title('IIR低通滤波器滤波后的频谱')
xlabel('频率/Hz');
ylabel('幅值');
[x1,fs,bits]=wavread('录音.wav'); %读取语音信号的数据，赋给变量x1 y1=fft(x1,1024); %对采样后信号做1024点FFT变换 f=fs*(0:511)/1024;
figure(1)
subplot(2,1,1);
plot(f,abs(y1(1:512)));
title('原始语音信号频谱')
xlabel('频率/Hz');
ylabel('幅值')
subplot(2,1,2);
plot(abs(y1(1:1024))) %采样后信号的FFT频谱图
title('原始语音信号FFT频谱')
xlabel('点数N');
ylabel('幅值');
[x1,fs,bits]=wavread('录音.wav'); %读取语音信号的数据，赋给变量x1 figure(1)
plot(x1) %做原始语音信号以44.1k采样后的时域图形 title('原始语音采样后时域信号');
xlabel('时间轴 n');
ylabel('幅值 A');
figure(2)
freqz(x1) %绘制原始语音信号采样后的频率响应图 title('原始语音信号采样后频率响应图');
[x1,fs,bits]=wavread('录音.wav'); %读取语音信号的数据，赋给变量x1 sound(x1,fs); %播放语音信号