语音信号处理与滤波(数字信号处理课程设计)

《数字信号处理》

课程设计报告

题 目： 语音信号处理与滤波 专 业： 通信工程 学 号： 姓 名： 指导教师： *** 成 绩：

二〇一二 年 一 月 三 日

一、课程设计目的

综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解。

二、课程设计基本要求及分组 1、对知识点的掌握要求：

（1）熟悉离散信号和系统的时域特性。 （2）掌握序列快速傅里叶变换FFT方法。

（3）学会MATLAB的使用，掌握MATLAB的程序设计方法。 （4）利用MATLAB对语音信号进行频谱分析。 （5）掌握MATLAB设计FIR和IIR数字滤波器的方法。 三、详细设计过程

1、 问题描述：录制一段自己的语音信号，取不同的数据点对语音信号进行频谱

分析；对所有数据进行插值和抽取处理，改变抽样率再对信号进行频谱分析；设计FIR和IIR数字滤波器，并对被抽样后的语音信号进行滤波，分析滤波后信号的时域和频域特征，回放语音信号。

2、详细操作步骤和部分运行结果 （1）采集语音信号：

利用Windows下的录音机，录制一段自己的话音，时间控制在4s左右；然

后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过使用wavread函数，理解采样频率、采样位数等概念。

wavread函数调用格式如下：

y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。 [y,fs,nbits]=wavread(file)，采样值放在向量y中，fs表示采样频率(Hz)，

bits表示采样位数。

y=wavread(file,N)，读取前N点的采样值放在向量y中。

y=wavread(file,[N1，N2])，读取从N1点到N2点的采样值放在向量y中。 录制并存储语音的程序代码：(函数名'luyin'被主函数‘main’调用)

fs=8000; t=4;

a=msgbox('请说话...');

x=wavrecord(t*fs,fs,1,'Double');%录制语音信号 close(a);

wavwrite(x,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计.wav');% 存储录音信号

原始语音时域信号波形图：

程序代码：

[x1,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 \\课程设计.wav'); sound(x1,fs,bits);%播放语音 figure(1);

subplot(1,1,1);plot(x1); %做原始语音信号的时域图形 title('原始语音信号 作品');

xlabel('时间 n');ylabel('音量 n'); 运行结果：

（2）对语音信号进行频谱分析：

在MATLAB中，利用函数fft对信号进行快速傅里叶变换，得到信号的频

谱特性。首先画出语音信号的时域波形，然后对语音信号进行频谱分析。

程序代码：

[x1,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计

忠\\课程设计.wav'); y1=fft(x1,1024); f=fs*(0:511)/1024; figure(1)

y1=fft(x1); %做length(x1)点的FFT

y1=fftshift(y1);%平移，是频率中心为0

derta_fs = fs/length(x1);%设置频谱的间隔，分辨率

subplot(2,1,1);plot([-fs/2:derta_fs: fs/2-derta_fs],abs(y1)); %画出原始语音信号的频谱图

title('原始语音信号的频谱 杜作品');

subplot(2,1,2); plot(abs(y1(1:1024))); title('原始语音信号FFT频谱 杜作品'); xlabel('点数N'); ylabel('幅值

原始语音信号频谱和FFT转换后的频谱图

（3）设计数字滤波器和画出频率响应：

根据语音信号的特点给出有关滤波器的性能指标：

① 低通滤波器性能指标：fp=1000Hz，fc=1200Hz，As=100dB，Ap=1dB； ② 高通滤波器性能指标：fc=4800Hz，fp=5000Hz，As=100dB，Ap=1dB； ③ 带通滤波器性能指标：fp1=1200Hz，fp2=3000Hz，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB。

首先用窗函数法（矩形窗(Rectangular window)、三角窗(Triangular window)、

汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。）设计上面要求的三种滤波器。在MATLAB中，利用函数fir1设计FIR滤波器；然后再用双线性变换法设计上面要求的三种滤波器，利用函数butte、cheby1和ellip设计IIR滤波器；最后，利用MATLAB中的函数freqz画出各种滤波器的频率响应。

（4）IIR滤波器的设计过程

IIR低通滤波器： 程序代码：()

fs=8000;

T=1/fs;

duration=3;

[y,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设 计.wav'); N=length(y); t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1); plot(y);

title('原始信号时域的图形') xlabel('时间/t') ylabel('幅值') subplot(2,1,2); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形') xlabel('频率/Hz')

ylabel('幅值/db')

fs=8000; fp=1000; fst=1200; ap=1; ast=100; wp=2*fp/fs; wst=2*fst/fs;

[N,wc]=cheb1ord(wp,wst,ap,ast); [num,den]=cheby1(N,ap,wc,'low'); [h,w]=freqz(num,den); subplot(2,1,1);

plot(w/pi,abs(h));grid; xlabel('\\omega/\\pi');

ylabel('振幅（幅值）');

title('契比雪夫Ⅰ型低通滤波器的幅频响应'); subplot(2,1,2);

plot(w/pi,20*log10(abs(h)));grid; xlabel('\\omega/\\pi');

ylabel('振幅(分贝)');

title('契比雪夫Ⅰ型低通滤波器的幅频响应');

y0=filter(num,den,y); fftwav=fft(y0); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y0);

title('通过契比雪夫Ⅰ型低通滤波器后的信号图形') xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz') ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1))));

title('滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y0,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 \\ 课程设计IIR_LP.wav');

[y0,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计IIR_LP.wav');%语音回放 soundview(y0,fs,'课程设计IIR_LP.wav

运行结果如下图：

IIR高通滤波器设计

程序代码：

fs=8000; T=1/fs;

duration=3;

[y,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜艺忠\\课程设 计.wav'); N=length(y); t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1); plot(y);

title('原始信号时域的图形') xlabel('时间/t') ylabel('幅值') subplot(2,1,2);

q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形') xlabel('频率/Hz') ylabel('幅值/db')

Fs=16000; fp=4000; fst=3800; ap=1; ast=100; wp=2*fp/Fs; wst=2*fst/Fs;

[N,wc]=cheb1ord(wp,wst,ap,ast); [num,den]=cheby1(N,ap,wc,'high'); [h,w]=freqz(num,den); subplot(2,1,1);

plot(w/pi,abs(h));grid; xlabel('\\omega/\\pi');

ylabel('振幅（幅值）');

title('契比雪夫Ⅰ型高通滤波器的幅频响应'); subplot(2,1,2);

plot(w/pi,20*log10(abs(h)));grid; xlabel('\\omega/\\pi');

ylabel('振幅(分贝)');

title('契比雪夫Ⅰ型高通滤波器的幅频响应');

y0=filter(num,den,y); fftwav=fft(y0); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y0);

title('通过契比雪夫Ⅰ型高通滤波器后的信号图形') xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz')

ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1))));

title('滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y0,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 \\课程 设计IIR_HP.wav');

[y0,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计IIR_HP.wav');%语音回放 soundview(y0,fs,'课程设计IIR_HP.wav');

程序运行结果如下图：

IIR带通滤波器的设计 程序代码:

fs=8000; T=1/fs;

duration=3;

[y,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设 计.wav'); N=length(y);

t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1); plot(y);

title('原始信号时域的图形') xlabel('时间/t') ylabel('幅值') subplot(2,1,2); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形') xlabel('频率/Hz')

ylabel('幅值/db')

fs=8000; ast=100; ap=1;

wp=[1200 3000]*2/fs; wst=[1000 3200]*2/fs;

[N,wn]=cheb1ord(wp,wst,ap,ast); [num,den]=cheby1(N,ap,wn); [h,w]=freqz(num,den); subplot(2,1,1);

plot(w/pi,abs(h));grid; xlabel('\\omega/\\pi');

ylabel('振幅（幅值）');

title('契比雪夫Ⅰ型带通滤波器的幅频响应'); subplot(2,1,2);

plot(w/pi,20*log10(abs(h)));grid; xlabel('\\omega/\\pi');

ylabel('振幅(分贝)');

title('契比雪夫Ⅰ型带通滤波器的幅频响应'); y1=filter(num,den,y); fftwav=fft(y1); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y1);

title('通过IIR带通滤波器后的信号图形')

xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz') ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1))));

title('带通IIR滤波器滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y1,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程 设计IIR_BP.wav');

[y1,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计IIR_BP.wav');%语音回放) soundview(y1,fs,'课程设计IIR_BP.wav');

程序运行结果：

（5）FIR滤波器的设计

FIR低通滤波器：

程序代码：

fs=8000; T=1/fs;

duration=3;

[y,fs,sbits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜艺\\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设 计.wav'); N=length(y); t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1); plot(y);

title('原始信号时域的图形 ') xlabel('时间/t') ylabel('幅值 ') subplot(2,1,2); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形') xlabel('频率/Hz') ylabel('幅值/db')

fs=8000; fp=1000; fb=1200; ap=1; ast=100;

wp=2*pi*fp/fs; wst=2*pi*fb/fs; deltaW=wst-wp; wc=(wst+wp)/2;

N=ceil(2*6.4*pi/deltaW);

hdWindow=ideallp(wc,N);%---ideallp（）函数（非系统自有函数）在系 统安装目录的WORK子目录中，ideallp.m wdWindow=Kaiser(N,10.056); hr=wdWindow.*hdWindow'; n=0:N-1;

subplot(2,2,1);

stem(n,wdWindow); xlabel('时间')

ylabel('振幅 作品') title('凯泽窗') subplot(2,2,2); stem(n,hr);

xlabel('时间') ylabel('振幅')

title('凯泽窗冲激响应') [H,W]=freqz(hr,1); subplot(2,2,3); plot(W/pi,abs(H)); xlabel('\\omega/\\pi')

ylabel('绝对振幅')

title('低通滤波器幅频特性（1）') subplot(2,2,4);

plot(W/pi,20*log10(abs(H))); xlabel('\\omega/\\pi')

ylabel('相对振幅')

title('低通滤波器幅频特性（2）')

y0=fftfilt(hr,y); fftwav=fft(y0); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y0);

title('通过凯泽低通滤波器后的信号图形') xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz') ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1)))); title('滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y0,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 艺\\ 课程设计FIR_LP.wav');

[y0,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计FIR_LP.wav');%语音回放 soundview(y0,fs,'课程设计FIR_LP.wav');

程序运行结果：

FIR高通滤波器：

程序代码：

fs=8000; T=1/fs;

duration=3;

[y,fs,sbits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计.wav');

N=length(y); t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1);

plot(y);

title('原始信号时域的图形') xlabel('时间/t') ylabel('幅值') subplot(2,1,2); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形') xlabel('频率/Hz') ylabel('幅值/db')

Fs=16000; fp=4000; fst=3800; ast=100; ap=1;

wp=2*pi*fp/Fs; wst=2*pi*fst/Fs; wc=(wst+wp)/2; deltaW=wp-wst;

N0=ceil(2*6.4*pi/deltaW); %---查表7-3（P342）凯泽窗 N=N0+mod(N0+1,2); %---确保N为奇数 N=ceil(2*6.4*pi/deltaW);

hdWindow=ideallp(pi,N)-ideallp(wc,N);%---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，ideallp.m

wdWindow=Kaiser(N,10.056); hr=wdWindow.*hdWindow'; n=0:N-1;

subplot(2,2,1); stem(n,wdWindow);

xlabel('时间') ylabel('振幅') title('凯泽窗') subplot(2,2,2); stem(n,hr);

xlabel('时间') ylabel('振幅')

title('凯泽窗冲激响应') [H,W]=freqz(hr,1); subplot(2,2,3); plot(W/pi,abs(H)); xlabel('\\omega/\\pi')

ylabel('绝对振幅')

title('凯泽高通滤波器幅频特性（1）')

subplot(2,2,4);

plot(W/pi,20*log10(abs(H))); xlabel('\\omega/\\pi')

ylabel('相对振幅')

title('凯泽高通滤波器幅频特性（2）')

y0=fftfilt(hr,y); fftwav=fft(y0); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y0);

title('通过凯泽高通滤波器后的信号图形') xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz') ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1))));

title('凯泽高通滤波器滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y0,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程 设计FIR_HP.wav');

[y0,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 忠\\课程设计FIR_HP.wav');%语音回放

soundview(y0,fs,'课程设计FIR_HP.wav');

程序运行结果：

FIR带通滤波器：

程序代码：

fs=8000; T=1/fs;

duration=3;

[y,fs,sbits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 \\课程设计.wav');

soundview('D:\\MATLAB\\R2010b\\work\\数字信号课程设计 杜\\课程设计.wav');

N=length(y); t=(1:N)/fs; df=fs/N; n1=1:N/2; f=(n1-1)*df; subplot(2,1,1); plot(y);

title('原始信号时域的图形') xlabel('时间/t') ylabel('幅值') subplot(2,1,2); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('原始信号频谱图形')

xlabel('频率/Hz') ylabel('幅值/db')

fs=8000; fp1=1200; fp2=3000; fst1=1000; fst2=3200; ast=100; ap=1;

wp1=2*pi*fp1/fs; wp2=2*pi*fp2/fs; wst1=2*pi*fst1/fs; wst2=2*pi*fst2/fs; wc1=(wst1+wp1)/2; wc2=(wst2+wp2)/2;

deltaW=min((wp1-wst1),(wst2-wp2));

N0=ceil(2*6.4*pi/deltaW); %---查表7-3（P342）凯泽窗 N=N0+mod(N0+1,2); %---确保N为奇数 N=ceil(2*6.4*pi/deltaW);

hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%---ideallp（）函数（非系统 自有函数）在系统安装目录的WORK子目录中，ideallp.m wdWindow=Kaiser(N,10.056); hr=wdWindow.*hdWindow'; n=0:N-1;

subplot(2,2,1); stem(n,wdWindow); xlabel('时间') ylabel('振幅') title('凯泽窗') subplot(2,2,2); stem(n,hr);

xlabel('时间') ylabel('振幅')

title('凯泽窗冲激响应') [H,W]=freqz(hr,1); subplot(2,2,3); plot(W/pi,abs(H)); xlabel('\\omega/\\pi')

ylabel('绝对振幅')

title('带通滤波器幅频特性（1）') subplot(2,2,4);

plot(W/pi,20*log10(abs(H))); xlabel('\\omega/\\pi')

ylabel('相对振幅')

title('带通滤波器幅频特性（2）')

y0=fftfilt(hr,y); fftwav=fft(y0); subplot(2,2,1); plot(y);

title('原始信号的波形') xlabel('时间/t') ylabel('幅值') subplot(2,2,2); plot(y0);

title('通过带通滤波器后的信号图形') xlabel('时间/t') ylabel('幅值') subplot(2,2,3); q=fft(y);

plot(f,20*log10(abs(q(n1)))); title('滤波前的信号频谱图') xlabel('频率/Hz') ylabel('幅值/db') subplot(2,2,4);

plot(f,20*log10(abs(fftwav(n1))));

title('带通滤波器滤波后信号的频谱图') xlabel('频率/Hz') ylabel('幅值/db')

wavwrite(y0,fs,'D:\\MATLAB\\R2010b\\work\\数字信号课程设计 \\课程 设计FIR_BP.wav');

[y0,fs,bits]=wavread('D:\\MATLAB\\R2010b\\work\\数字信号课程设计

忠\\课程设计FIR_BP.wav');%语音回放 soundview(y0,fs,'课程设计FIR_BP.wav');

程序运行结果：

杜艺

（6）设计系统界面

在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤

波器的频率响应，选择信号等。 主函数程序（main）：

b=menu('忠作品 版权所有，翻版必究','录音','原始信号采样后时域图和频谱 图','FIR滤波器','IIR滤波器','退出'); if b==5 b=0; End if b==1 luyin; main;

else if b==2

temp=menu('作品 版权所有，翻版必究','播放原始语音并生成时域图',' 原始语音频谱图及做FFT变换后频谱图','返回'); if temp==1 cjyy; main;

else if temp==2 cppfx; main; else

main; end end

else if b==3

temp=menu('忠作品 版权所有，翻版必究','FIR低通滤波器','FIR高 通滤波器','FIR带通滤波器','返回'); if temp==1 FIR_LP; main; else if temp==2 FIR_HP; main; else if temp==3 FIR_BP; main; else main; end end end else if b==4

temp=menu('杜作品 版权所有，翻版必究','IIR低通滤波器','IIR高 通滤波器','IIR带通滤波器','返回'); if temp==1 IIR_LP; main;

else if temp==2 IIR_HP; main;

else if temp==3 IIR_BP; main; else

main; end end end end end end end

界面如下图所示：

四、调试分析

1、在调试过程中发现原先安装的MATLAB7.0有问题，缺少切比雪夫函数，后来重装

MATLAB2010版的，运行IIR滤波器程序顺利。

2.、在MATAB中没有ideallp（）函数（非系统自有函数），在网上下载个放在系统

五、心得体会

通过一学期的数字信号处理的学习，在这次做课程设计的过程中，真的发现自己没学到啥东西，所以导致做的过程中困难重重，拿着课本只好重头看，通过网上查资料，同学之间的相互合作，最终我们还是战胜了困难，顺利完成了此

安装目录的work子目录中，FIR滤波器程序能运行。

次的课程设计。我深深的体会到了知识的重要性，不管是为了应付考试还是怎么样也好，对于自身的帮助还是很大的，整整两天坐在电脑面前的功夫没有白费，我还是收获了很多，对于Matlab软件的使用也更加熟悉了。

我觉得这才是真正学到知识的，也锻炼了自己的动手实践能力，耐心和毅力。学会团结，才能把事情做得更好。

七、参考文献和附录

[1]丁美玉，高西全 数字信号处理2版。西安：西安电子科技大学出版社，2001 [2]怀琛 数字信号处理教程——MATLAB释疑与实现。北京：电子工业出版社，2004 [3]王宏 MATLAB6.5及其在信号处理中的应用。北京：清华大学出版社，2004 [4]刘顺兰，吴杰 数字信号处理。西安：西安电子科技大学出版社，2003