数字信号处理课程设计..

课程设计报告
课程名称：数字信号处理
课题名称：语音信号的处理与滤波姓名：
学号：
院系：
专业班级：
指导教师：
完成日期： 2013年7月2日
目录
第1部分课程设计报告 (3)
一．设计目的 (3)
二．设计内容 (3)
三．设计原理 (3)
四．具体实现 (5)
1.录制一段声音 (5)
2.巴特沃斯滤波器的设计 (8)
3.将声音信号送入滤波器滤波 (13)
4.语音信号的回放 (19)
5.男女语音信号的频谱分析 (19)
6.噪声的叠加和滤除 (22)
五．结果分析 (27)
第2部分课程设计总结 (28)
一．参考文献 (28)
第1部分课程设计报告
一．设计目的
综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

二．设计内容
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除；
三．设计原理
1.在Matlab软件平台下，利用函数wavrecord()，wavwrite(),wavread(),wavplay()对语音信号进行录制，存储，读取，回放。

2.用y=fft(x)对采集的信号做快速傅立叶变换，并用
[h1,w]=freqz(h)进行DTFT 变换。

3.掌握FIR DF 线性相位的概念，即线性相位对)()(ωH n h 、及零
点的约束，了解四种FIR DF 的频响特点。

4.在Matlab 中，FIR 滤波器利用函数fftfilt 对信号进行滤
波。

5.抽样定理
连续信号经理想抽样后时域、频域发生的变化（理想抽样信号
与连续信号频谱之间的关系）
理想抽样信号能否代表原始信号、如何不失真地还原信号即由
离散信号恢复连续信号的条件（抽样定理）
理想采样过程描述：
时域描述：
频域描述：利用傅氏变换的性质，时域相乘频域卷积，若
()()T T j t δ∆Ω↔ 则有
1ˆ()()()2a a T X j X j j πΩ=Ω*∆Ω ˆ()a
X j Ω与()a X j Ω的关系：理想抽样信号的频谱是连续信号频谱的ˆ()()()()()()()a a T a a n n x t x t t x t t nT x nT t nT δδδ∞∞=-∞=-∞==-=-∑∑()()T n t t nT δδ∞=-∞=
-∑ˆ()()a a X j x t Ω↔()()
a a X j x t Ω↔121ˆ()()()a a a s k k X j X j jk X j jk T T T π∞∞=-∞=-∞
Ω=Ω-=Ω-Ω∑∑
周期延拓,重复周期为Ωs(采样角频率)。

如果：
即连续信号是带限的，且信号最高频率不超过抽样频率的二分之一，则可不失真恢复。

奈奎斯特采样定理：要使实信号采样后能够不失真还原，采样频率必须大于信号最高频率的两倍：2 2s h s h f f Ω≥Ω≥或
四．具体实现
1.录制一段声音
1.1录制并分析
在MATLAB 中用wavrecord 、wavread 、wavplay 、wavwrite 对声音进行录制、读取、回放、存储。

程序如下：
Fs=8000; %抽样频率
time=3; %录音时间
fprintf('按Enter 键录音%ds',time); %文字提示
pause; %暂停命令
fprintf('录音中......');
x=wavrecord(time*Fs,Fs,'double'); %录制语音信号 fprintf('录音结束'); %文字提示
fprintf('按Enter 键回放录音');
pause; %暂停命令
()/2
()0/2
a
s a s X j X j ⎧ΩΩ<Ω⎪Ω=⎨Ω≥Ω⎪⎩
wavplay(x,Fs); %按任意键播放语音信号
wavwrite(x,Fs,'C:\Users\acer\Desktop\数字信号
\sound.wav'); %存储语音信号
N=length(x);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
figure(2);
subplot(2,1,1);
plot(x);%录制信号的时域波形
title('原始信号的时域波形');%加标题
ylabel('幅值/A');%显示纵坐标的表示意义
grid;%加网格
y0=fft(x);%快速傅立叶变换
figure(2);
subplot(2,1,2);
plot(f,abs(y0(n1)));%原始信号的频谱图
title('原始信号的频谱图');%加标题
xlabel('频率w/pi');%显示横坐标表示的意义
ylabel('幅值 ');%显示纵坐标表示的意义
title('原始信号的频谱图');%加标题
grid;%加网格
图1.1 原始信号的时域与频谱图
1.2滤除无效点
针对实际发出声音落后录制动作半拍的现象，如何拔除对无效点的采样的问题: 出现这种现象的原因主要是录音开始时，人的反应慢了半拍，导致出现了一些无效点，而后而出现的无效的点，主要是已经没有声音的动作，先读取声音出来，将原始语音信号时域波形图画出来，根据己得到的信号，可以在第二次读取声音的后面设定采样点，取好有效点，画出滤除无效点后的语音信号时域波形图，对比可以看出。

这样就可以解决这个问题。

x=wavread('C:\Users\acer\Desktop\数字信号\sound.wav',
[4000,24000] );
%从4000点截取到24000结束
plot (x);%画出截取后的时域图形
title('截取后的声音时域图形');%标题
xlabel('频率');ylabel('振幅') ;
grid; %画网格
图1.2 去除无效点
2.巴特沃斯滤波器的设计
2.1设计巴特沃思低通滤波器
MATLAB程序如下。

滤波器图如图3.3所示。

%低通滤波
fp=1000;fs=1200;Fs=22050;
rp=1;rs=100;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Fs1=1;
wap=2*tan(wp/2);
was=2*tan(ws/2);
[N,wc]=buttord(wap,was,rp,rs,'s');
[B,A]=butter(N,wc,'s');
[Bz,Az]=bilinear(B,A,Fs1);
figure(1);
[h,w]=freqz(Bz,Az,512,Fs1*22050);
plot(w,abs(h));
title('巴特沃斯低通滤波器');
xlabel('频率（HZ）');ylabel('耗损（dB）'); gridon;
图2.1 巴特沃思低通滤波器
2.2设计巴特沃思高通滤波器
MATLAB程序如下。

滤波器图如图3.5所示。

%高通滤波
fp=4800;fs=5000;Fs=22050;
rp=1;rs=100;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
T=1;Fs1=1;
wap=2*tan(wp/2);
was=2*tan(ws/2);
[N,wc]=buttord(wap,was,rp,rs,'s');
[B,A]=butter(N,wc,'high','s');
[Bz,Az]=bilinear(B,A,Fs1);
figure(1);
[h,w]=freqz(Bz,Az,512,Fs1*22050);
plot(w,abs(h));
title('巴特沃斯高通滤波器');
xlabel('频率（HZ）');ylabel('耗损（dB）'); grid on;
图2.2巴特沃思高通滤波器
2.3设计巴特沃思带通滤波器
MATLAB程序如下。

滤波器图如图3.7所示。

%带通滤波
fp=[1200,3000];fs=[1000,3200];Fs=8000;
rp=1;rs=100;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
T=1;Fs1=1;
wap=2*tan(wp/2);
was=2*tan(ws/2);
[N,wc]=buttord(wap,was,rp,rs,'s');
[B,A]=butter(N,wc,'s');
[Bz,Az]=bilinear(B,A,Fs1);
figure(4);
[h,w]=freqz(Bz,Az,512,Fs1*1000);
plot(w,abs(h));
title('巴特沃斯带通滤波器');
xlabel('频率（HZ）');ylabel('耗损（dB）'); grid on;
图2.3巴特沃思带通滤波器
3．将声音信号送入滤波器滤波
x=wavread('C:\Users\acer\Desktop\数字信号\sound.wav');%播放原始信号
wavplay(x,fs); %播放原始信号
N=length(x);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
figure(4);
subplot(4,2,1);
plot(x);%录制信号的时域波形
title('原始信号的时域波形');%加标题ylabel('幅值/A');%显示纵坐标的表示意义grid;%加网格
y0=fft(x);%快速傅立叶变换
subplot(4,2,3);
plot(f,abs(y0(n1)));%原始信号的频谱图
title('原始信号的频谱图');%加标题
xlabel('频率w/pi');%显示横坐标表示的意义ylabel('幅值 ');%显示纵坐标表示的意义
title('原始信号的频谱图');%加标题
grid;%加网格
3．1低通滤波器滤波
fs=8000;
beta=10.056;
wc=2*pi*1000/fs;
ws=2*pi*1200/fs;
width=ws-wc;
wn=(ws+wc)/2;
n=ceil(12.8*pi /width);
h=fir1(n,wn/pi,'band',kaiser(n+1,beta)); [h1,w]=freqz(h);
ys=fftfilt(h,x);%信号送入滤波器滤波，ys为输出
fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换figure(4);
subplot(4,2,2);%在四行两列的第二个窗口显示图形
plot(ys);%信号的时域波形
title('低通滤波后信号的时域波形');%加标题
xlabel('频率w/pi');
ylabel('幅值/A');%显示标表示的意义
grid;%网格
subplot(4,2,4); %在四行两列的第四个窗口显示图形
plot(f, abs(fftwave(n1)));%绘制模值
xlabel('频率w/pi');
ylabel('幅值/A');%显示标表示的意义
title('低通滤波器滤波后信号的频谱图');%标题
grid;%加网格
wavplay(ys,8000); %播放滤波后信号
3．2高通滤波器滤波
fs=8000;
beta=10.056;
ws=2*5000/fs;
wc=2*4800/fs;
width=ws-wc;
wn=(ws+wc)/2;
n=ceil(12.8*pi/width);
h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));
[h1,w]=freqz(h);
ys=fftfilt(h,x);%将信号送入高通滤波器滤波
subplot(4,2,5);%在四行两列的第五个窗口显示图形
plot(ys);%信号的时域波形
xlabel('频率w/pi');
ylabel('幅值/A');%显示标表示的意义
title('高通滤波后信号的时域波形');%标题
ylabel('幅值/A');%显示纵坐标的表示意义
grid;%网格
fftwave=fft(ys); %将滤波后的语音信号进行快速傅立叶变换subplot(4,2,7);%在四行两列的第七个窗口显示图形
plot(f,abs(fftwave(n1)));%绘制模值
axis([0 1 0 50]);
xlabel('频率w/pi');
ylabel('幅值/A');%显示标表示的意义
title('高通滤波器滤波后信号的频谱图');%标题
grid;%加网格
wavplay(ys,8000); %播放滤波后信号
3．3带通滤波器
fs=8000;
beta=10.056;
wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;
ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;
width=ws1-wc1;
wn1=(ws1+wc1)/2; wn2=(ws2+wc2)/2;
wn=[wn1 wn2];
n=ceil(12.8/width*pi);
h=fir1(n,wn/pi,'band',kaiser(n+1,beta));
[h1,w]=freqz(h);
ys1= fftfilt(h,x); %将信号送入高通滤波器滤波
figure(4);
subplot(4,2,6);%在四行两列的第六个窗口显示图形
plot(ys1);%绘制后信号的时域的图形
title('带通滤波后信号的时域波形');%加标题
xlabel('频率w/pi');
ylabel('幅值/A');%显示纵坐标表示的意义
grid;%网格
fftwave=fft(ys1);%对滤波后的信号进行快速傅立叶变换subplot(4,2,8);%在四行两列的第八个窗口显示图形
plot(f, abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);
xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('带通滤波器滤波后信号的频谱图');%加标题 grid;%网格
wavplay(ys1,8000); %播放滤波后信号 图形如下：
1
2
3
x 10
4
-10
1
原始信号的时域波形
幅值/A
0.510
100
200原始信号的频谱图
频率w/pi
幅值
123
x 10
4
-0.50
0.5低通滤波后信号的时域波形频率w/pi
幅值/A
00.51
0100
200频率w/pi
幅值/A
低通滤波器滤波后信号的频谱图0
123
x 10
4
-0.2
0.2频率w/pi
幅值/A
高通滤波后信号的时域波形
0.51
50
频率w/pi
幅值/A
高通滤波器滤波后信号的频谱图0123
x 10
4-0.50
0.5带通滤波后信号的时域波形频率w/pi
幅值/A
0.51
50
频率w/pi
幅值/A
带通滤波器滤波后信号的频谱图
分析：三个滤波器滤波后的声音与原来的声音都发生了变化。

其中低
通的滤波后与原来声音没有很大的变化，其它两个都又明显的变化
4.语音信号的回放
sound(xlow,Fs,bits); %在Matlab中，函数sound可以对声音进行回放，其调用格式:
sound(xhigh, Fs,bits); %sound (x, Fs, bits);
sound(xdaitong, Fs,bits);
5.男女语音信号的频谱分析
5.1 录制一段异性的声音进行频谱分析
Fs=8000; %抽样频率
time=3; %录音时间
fprintf('按Enter键录音%ds',time); %文字提示
pause; %暂停命令
fprintf('录音中......');
x=wavrecord(time*Fs,Fs,'double'); %录制语音信号fprintf('录音结束'); %文字提示
fprintf('按Enter键回放录音');
pause; %暂停命令
wavplay(x,Fs); %按任意键播放语音信号
wavwrite(x,Fs,'C:\Users\acer\Desktop\数字信号
\sound2.wav'); %存储语音信号
5.2 分析男女声音的频谱
x=wavread(' C:\Users\acer\Desktop\数字信号\sound2.wav ');%播放原始信号,解决落后半拍
wavplay(x,fs); %播放原始信号
N=length(x);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
figure(1);
subplot(2,2,1);
plot(x);%录制信号的时域波形
title('原始女生信号的时域波形');%加标题
ylabel('幅值/A');%显示纵坐标的表示意义
grid;%加网格
y0=fft(x);%快速傅立叶变换
subplot(2,2,2);
plot(f,abs(y0(n1)));%原始信号的频谱图
title('原始女生信号的频谱图');%加标题
xlabel('频率w/pi');%显示横坐标表示的意义
ylabel('幅值 ');%显示纵坐标表示的意义
grid;%加网格
[y,fs,bits]=wavread(' C:\Users\acer\Desktop\数字信号\sound.wav ');% 对语音信号进行采样
wavplay(y,fs); %播放原始信号
N=length(y);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
subplot(2,2,3);
plot(y);%录制信号的时域波形
title('原始男生信号的时域波形');%加标题
ylabel('幅值/A');%显示纵坐标的表示意义
grid;%加网格
y0=fft(y);%快速傅立叶变换
subplot(2,2,4);%在四行两列的第三个窗口显示图形
plot(f,abs(y0(n1)));%原始信号的频谱图
title('原始男生信号的频谱图');%加标题
xlabel('频率w/pi');%显示横坐标表示的意义
ylabel('幅值 ');%显示纵坐标表示的意义
grid;%加网格
5.3男女声音的频谱图
1
2
3
x 10
4
-1
-0.50
0.5原始女生信号的时域波形幅值/A
0.51
50100
150原始女生信号的频谱图
频率w/pi
幅值
1
2
3
x 10
4
-0.50
0.5
原始男生信号的时域波形幅值/A
0.51
100
200
300原始男生信号的频谱图
频率w/pi
幅值
图5.3男女声音信号波形与频谱对比
分析：就时域图看，男生的时域图中振幅比女生的高，对于频谱图女生的高频成分比较多
6.噪声的叠加和滤除
6.1录制一段背景噪声
Fs=8000; %抽样频率 time=3; %录音时间 fprintf('按Enter 键录音%ds',time); %文字提示 pause; %暂停命令 fprintf('录音中......');
x=wavrecord(time*Fs,Fs,'double'); %录制语音信号
fprintf('录音结束'); %文字提示
fprintf('按Enter键回放录音');
pause; %暂停命令
wavplay(x,Fs); %按任意键播放语音信号wavwrite(x,Fs,'C:\Users\acer\Desktop\数字信号\噪
音.wav'); %存储语音信号
6.2 对噪声进行频谱的分析
[x1,fs,bits]=wavread(' C:\Users\acer\Desktop\数字信号\噪音.wav ');%对语音信号进行采样
wavplay(x1,fs); %播放噪声信号
N=length(x1);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
figure(5);
subplot(3,2,1);
plot(x1);%信号的时域波形
title('噪声信号的时域波形');
grid;
ylabel('幅值/A');
y0=fft(x1);%快速傅立叶变换
subplot(3,2,2);
plot(f,abs(y0(n1)));%噪声信号的频谱图
ylabel('幅值');
title('噪声信号的频谱图');
6.3原始信号与噪音的叠加
fs=8000;
[x,fs,bits]=wavread(' C:\Users\acer\Desktop\数字信号\sound.wav '); %对录入信号进行采样
[x1,fs,bits]=wavread(' C:\Users\acer\Desktop\数字信号\噪音.wav '); %对噪声信号进行采样
yy=x+x1;%将两个声音叠加
6.4叠加信号的频谱分析：
wavplay(yy,fs); %播放叠加后信号
N=length(yy);%返回采样点数
df=fs/N;%采样间隔
n1=1:N/2;
f=[(n1-1)*(2*pi/N)]/pi;%频带宽度
figure(5);
subplot(3,2,3);
plot(yy,'LineWidth',2);%信号的时域波形
title('叠加信号的时域波形');
xlabel('时间/t');
ylabel('幅值/A');
grid;
y0=fft(yy);%快速傅立叶变换
subplot(3,2,4);
plot(f,abs(y0(n1)));%叠加信号的频谱图
title('叠加信号的频谱图');
xlabel('频率w/pi');
ylabel('幅值/db');
grid;
6.5 设计一个合适的滤波器将噪声滤除
fs=18000;%采样频率
Wp=2*1000/fs;%通带截至频率
Ws=2*2000/fs;%阻带截至频率
Rp=1;%最大衰减
Rs=100;%最小衰减
[N,Wn]=buttord(Wp,Ws,Rp,Rs); %buttord函数（n为阶数，Wn 为截至频率）
[num,den]=butter(N,Wn); %butter函数（num为分子系数den 为分母系数）
[h,w]=freqz(num,den);%DTFT变换
ys=filter(num,den,yy);%信号送入滤波器滤波，ys为输出fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换figure(5);
subplot(3,2,5);
plot(ys);%信号的时域波形
title('低通滤波后信号的时域波形');%加标题
ylabel('幅值/A');%显示标表示的意义
grid;%网格
subplot(3,2,6);
plot(f, abs(fftwave(n1)));%绘制模值
title('低通滤波器滤波后信号的频谱图');%标题
xlabel('频率w/pi');
ylabel('幅值/A');%显示标表示的意义
grid;%加网格
wavplay(ys,8000); %播放滤波后信号
grid;
图形如下：
1
2
3
x 10
4
-1
1
噪声信号的时域波形幅值/A
00.5
1
50
100
幅值
噪声信号的频谱图
12
3
x 10
4
-10
1叠加信号的时域波形时间/t
幅值/A
0.51
0100
200
叠加信号的频谱图
频率w/pi
幅值/d b 0
1
2
3
x 10
4
-0.5
0.5
低通滤波后信号的时域波形幅值/A
0.51
100
200
低通滤波器滤波后信号的频谱图频率w/pi
幅值/A
图6.1噪音的叠加与滤除前后频谱对比
7.结果分析
1.录制刚开始时，常会出现实际发出声音落后录制动作半拍，可在[x,fs,bits]=wavread('d:\matlav\work\womamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:\matlav\work\womamaaiwo.wav'，[100 10000])，窗的长度可根据需要定义。

2. 语音信号通过低通滤波器后，把高频滤除，声音变得比较低沉。

当通过高通滤波器后，把低频滤除，声音变得比较就尖锐。

通过带通滤波器后，声音比较适中。

3.通过观察男生和女生图像知：时域图的振幅大小与性别无关，只与说话人音量大小有关，音量越大，振幅越大。

频率图中，女生高
频成分较多。

4.叠加噪声后，噪声与原信号明显区分,但通过低通滤波器后，噪声没有滤除，信号产生失真。

原因可能为噪声与信号频率相近无法滤除。

第2部分课程设计总结通过本次课程设计，使我们对数字信号处理相关知识有了更深刻的理解，尤其是对各种滤波器的设计。

在设计的过程中遇到了很多问题，刚刚开始时曾天真的认为只要把以前的程序改了参数就可以用了，可是问题没有我想象中的那么简单，单纯的搬程序是不能解决问题的。

通过查阅资料和请教同学收获了很多以前不懂的理论知识。

再利用所学的操作，发现所写的程序还是没有能够运行，通过不断地调试，运行，最终得出了需要的结果。

整个过程中学到了很多新的知识，特别是对Matlab的使用终于有些了解。

在以后的学习中还需要深入了解这方面的内容。

在这次的课程设计中让我体会最深的是：知识来不得半点的马虎。

也认识到自己的不足，以后要进一步学习。

八．参考文献
[1]数字信号处理教程(第三版) 程佩青清华大学出版社
[2]MATLAB信号处理刘波文忠电子工业出版社
[3]MATLAB7.1及其在信号处理中的应用王宏
清华大学出版社
[4]MATLAB基础与编程入门张威西安电子科技大学出版社
[5] 数字信号处理及其MATLAB实验赵红怡张常
化学工业出版社
[6]MATLAB信号处理详解陈亚勇等人民邮电出版社
[7] 数字信号处理钱同惠机械工业出版社。