数字信号处理
课程设计报告
课题名称：语音信号的处理与滤波
系 别：
学 号：
姓 名：
班 级：
指导教师：
20##年10月25日
目录
摘录…………………………………………………………….3
正文…………………………………………………………….4
一、设计内容…………………………………………………..4
二、设计原理…………………………………………………..4
三、设计过程………………………………………………..…7
1、语音信号的录制…………………………………………....7
2、取点频谱分析……………………………………………...8
3、低通滤波器设计…………………………………………..10
4、滤滤前后信号对比………………………………………..12
结论…………………………………………………………...14
参考文献……………………………………………………...15
摘要
随着计算机和信息科学的飞速发展，信号处理逐渐发展成为一门独立的学科，成为信息科学的重要组成部分，在语音处理、雷达、图像处理、通信、生物医学工程等众多领域中得到广泛应用。 Matlab语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，Matlab功能强大、简单易学、编程效率高，深受广大科技工作者的喜爱。特别是Matlab还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。 数字信号处理课程在现代科学中具有很大重要性及自身特点，理解与掌握课程中的基本概念、基本原理、基本分析方法，对用Matlab进行数字信号处理课程设计的思路，具有很大帮助。 语音信号的处理与滤波是数字信号处理课程中常出现的课题，也是现代科学中值得深入究研的一个课题。 语音信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。
关键词：抽样频率；滤波器；抽样点；频谱分析
正文
一、设计内容：
1、熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数。
2．在MATLAB环境中，使用声音相关函数录制2秒左右自己的声音，抽样率是8000Hz/s。（考虑如何解决一个实际问题：录制刚开始时，常会出现实际发出声音落后录制动作半拍的现象，如何排除对这些无效点的采样？）
3．分别取8000个和16000个数据进行频谱分析，得到幅度和相位谱，比较二者异同并分析原因。
4．针对电话信道（最高3500Hz），设计一个FIR或IIR滤波器进行滤波，把抽样率转变为7000Hz/s，并进行频谱分析，得到幅度和相位谱。
5．把处理后的所有数据储存为声音文件，与原始声音进行比较。
二、设计原理
1、 语音信号的采样
l 语音信号的录制
语音信号的录制可通过MATLAB工具也可以通过windows自带的录音机录制一段自己的声音。在MATLAB里面，用函数wavrecord录制一段时间自定的声音，用函数wavwrite存储原始语音信号。其调用格式如下：
X=wavrecord(n ,fs ,ch ,dtype) ；
式中n 为抽样点数，决定了录音长度；
fs为抽样频率，默认值为11025Hz，还可根据要求自己选择合适的抽样率；
ch为声道数，默认值为1，表示单声道，如果指定为2，则抽样为双声道立体声数据；
dtype为抽样数据的存储格式，用字符串指定，可以是‘double’、‘single’、‘int16’、‘int8’，指定存储格式的同时也就规定了每个抽样值量化的精度，int8 对应8 位精度抽样，其它都是16 位抽样精度。
语音信号的存储可用函数wavwrite实现,其调用格式如下:
wavwrite(X,Fs,'f:\课程设计\sound.wav');
式中X为所录制的语音信号，Fs为抽样频率，引号内表示存储路径
l 语音信号的频谱分析
对语音信号进行频谱分析，是认识语音信号和处理语音信号的重要方法。下面将采用频谱来分析语音信号。
对一个语音时域离散信号{x(n)}，其频谱函数
是x(n)的傅里叶变换。傅里叶变换定义为：
其中j为虚数单位，
称为数字频率。
由于处理的信号是有限长的，即
，故实际采用的是离散傅里叶变换DFT 。 长度为N的序列x(n)，其DFT定义为：
，其中
如果x(n)为一个周期序列，得到的X(k)为x(n)的周期频谱。若x(n)不是周期序列，先对x(n)进行周期性扩展，把它看成某个周期性信号的一个周期，得到的X(k)是x(n) 频谱在一个周期的采样。
DFT有一种快速算法FFT，称为快速傅里叶变换。MATLAB中采用的就是FFT算法。
2、滤波器设计
滤波器的设计一般可分为无限冲激响应数字滤波器（IIR）和有限冲激响应数字滤波器（FIR）。本课题设计主要用FIR滤波器来设计。
l IIR数字滤波器
IIR数字滤波器的系统函数是Z的有理函数，可表示为系统函数的设计就是要确定系数或者零、极点，以使滤波器满足给定的性能要求。这种设计方法如下：
（1）用模拟滤波器理论来设计数字滤波器。在IIR数字滤波器的设计中较多采用这种方法。
（2）用计算机辅助设计，优化技术设计。即要求所设计的数字滤波器的输入x(n)等于x(t)的抽样信号x(Tn)时，输出y(n)也恰好等于y(t)的抽样信号。
l FIR数字滤波器
（1）FIR数字滤波器的算法描述
FIR数字滤波器又称卷积滤波器，因为它在时域上的输入x(n)和输出信号y(n)的数学关系是卷积运算：y(n)= hd(n)*x(n)。其中hd(n)是FIR滤波器的单位冲激响应函数，其表达式是由理想滤波器频响的幅度函数H（w）进行傅里叶逆变换得到的，式中wc为滤波器的截止频率。
（2）FIR滤波器的窗函数设计法
窗函数法的基本思想是，由于滤波器的时间响应函数hd(n)是无限长的，不能用FIR滤波的方法实现，所以要截取hd(n)幅度较大的部分，舍弃hd(n)幅度较小的部分来近似表示hd(n)。截取方法是用某一种窗函数和截取hd(n)的一段进行卷积以得到实际滤波器的相应hd(n)。窗函数不同，窗口宽度不同，实际频响会有较大区别。如加海明窗的频响曲线的通带和阻带特性比矩形窗好；窗口宽度N的阶数高，也能提高滤波器频响性能，但实时信号处理的时间会相应增加。因此，应根据滤波参数要求确定选用窗函数的类型和阶数。
（3）逼近
确定了技术指标后，就可以建立一个目标数字滤波器模型。通常采用理想的数字滤波器模型。之后，利用数字滤波器的设计方法，设计一个实际滤波器模型来逼近给定的目标。
（4）性能分析和计算机仿真
上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求，或利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。
（5）线性相位FIR滤波器具体的窗函数法设计方法
线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数w(n)，使其具有最窗宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特勒特(bartlett)窗、汉宁(hanning)窗、海明(hamming)窗、布莱克曼(blackman)窗和凯泽(kaiser)窗。
三、设计过程（程序实现）
l 语音信号的录制
l 语音信号的录制
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,'f:\课程设计\sound.wav'); %存储语音信号
在录音过程中，出现了这样一个现象：录制刚开始时，常会出现实际发出声音落后录制动作半拍的现象。出现这种现象的原因主要是录音开始时，人的反应慢了半拍，导致出现了一些无效点。要去掉这些无效点，可用函数wavread()来滤除。调用格式如下：
y= wavread('f:\课程设计\sound.wav',[3423 20050]);
subplot(2,1,1);
plot(y);
title('滤除无效点后的语音信号')
l 取8000点进行频谱分析
y=wavread('f:\课程设计\sound.wav'); %读取原始语音信号
u=fft(y,8000); %对语音信号做8000点的FFT变换
subplot(2, 2, 1); %分割窗口
plot(abs(u)); %绘出语音信号的幅度谱
axis([0, 4000 ,0,80]); %显示X轴、Y轴取值范围
title('8000点幅度谱'); %加图形说明
grid; %添加网格
subplot(2, 2, 2); %分割窗口
plot(angle(u)); %绘出语音信号的相位谱
axis([0, 4000 ,-4,4]); %显示X轴、Y轴取值范围
title('8000点相位谱'); %加图形说明
grid; %添加网格
l 取16000点进行频谱分析
y=wavread('f:\课程设计\sound.wav'); %读取原始语音信号
u=fft(y,16000); %对语音信号做16000点的FFT变换
subplot(2, 2, 3); %分割窗口
plot(abs(u)); %绘出语音信号的幅度谱
axis([0, 8000 ,0,150]); %显示X轴、Y轴取值范围
title('16000点幅度谱'); %加图形说明
grid; %添加网格
subplot(2, 2, 4); %分割窗口
plot(angle(u)); %绘出语音信号的相位谱
axis([0, 8000 ,-4,4]); %显示X轴、Y轴取值范围
title('16000点相位谱'); %加图形说明
grid %添加网格
8000点和16000点数的频谱分析图
取8000点和16000点进行频谱分析，可以发现两幅度谱和相位谱图形状大致相同，但在取数点的始段，16000点比8000点的振幅大，因为取16000点比取8000点数多，16000点显示得更详细，更接近真实值，而8000点相对于16000点来说只是大致趋向真实值而已。
l 滤波器的设计
function hd = ideallp(Wc,N); %自定义函数
% hd = ideallp(Wc,N)
% ----------------------------------------------
% hd = 点0 到 N-1之间的理想脉冲响应
% wc = 截止频率(弧度)
% N = 理想滤波器的长度
%
tao = (N-1)/2; % 理想脉冲响应的对称中心位置
n = [0: (N-1)]; % 设定脉冲响应长度
m = n - tao + eps; % 加一个小数以避免零作除数
hd = sin(Wc*m) ./ (pi*m); % 理想脉冲响应
function [db,mag,pha,W]=freqz_m(b,a); %自定义函数
%db为相对振幅(dB)；mag为绝对振幅；pha为相位响应；W频率样本点向量。
[H,W]=freqz(b,a,1000, 'whole');
%b和a分别为H(z)分子和分母多项式系数;
H=(H(1:501))';
W=(W(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
%。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
Wp=2* 3000/7000; Ap=1; %滤波器性能指标
Wr=2* 2400/7000; Ar=50;
Width= Wp-Wr; %过渡带宽度
Wc= (Wr+Wp)/2; %理想滤波器的截止频率
N1=ceil(6.6*pi/Width); %滤波器长度
N=N1+mod(N1+1,2); %确保N为奇数
n=0:1:N-1; %步长
hd=ideallp(Wc,N); %理想单位脉冲响应
window=hamming(N)'; %海明窗
h=window.*hd; %截取得到实际单位脉冲响应
[db,mag,pha,W]=freqz_m(h,1); %计算实际滤波器的幅频响应
subplot(2,1,1); %分割窗口
stem(n,h); %绘出实际单位脉冲响应
title('实际脉冲响应h(n)'); %加标题说明
grid; %添加网格
subplot(2,1,2); %分割窗口
plot(W/pi,db); %绘出实际滤波器的幅频响应
title('幅度响应(db)'); %加标题说明
grid %添加网格
%滤波后信号的频谱分析
fs=8000; %抽样频率
y=wavread('f:\课程设计\sound.wav'); %读取原始语音信号
y1=fftfilt(h,y); %将语音信号送入滤波器滤波，y1为输出
y2=fft(y1); %将滤波后的信号做FFT变换
subplot(2,1,1); %分割窗口
plot(abs(y2)); %绘出幅度图
axis([0, 12000 ,0,200]); %显示X轴、Y轴取值范围
title('滤波后信号的幅度谱图'); %加标题说明
grid; %添加网格
subplot(2,1,2); %分割窗口
plot(angle(y2)); %绘出相位图
axis([0, 12000 ,-4,4]); %显示X轴、Y轴取值范围
title('滤波后信号的相位谱图'); %加标题说明
grid; %添加网格
wavwrite(y1,fs,'f:\课程设计\sound1.wav'); %储存滤波后的语音信号
实际滤波器的幅度响应
滤波器的长度为121，最小阻带衰减53dB，满足滤波器设计要求。
滤波后信号的幅度、相位谱图
l 原始信号与滤波后的信号相比较
ys=wavread('f:\课程设计\sound.wav'); %读取原始语音信号
ys1= wavread('f:\课程设计\sound1.wav'); %读取滤波后语音信号
fprintf('按Enter键播放原始语音信号'); %文字提示
pause; %暂停命令
wavplay(ys,fs); %播放原始语音信号
fprintf('按Enter键播放滤波后语音信号'); %文字提示
pause; %暂停命令
wavplay(ys1,fs); %播放滤波后语音信号
yf=fft(ys); %对原始声音进行FFT变换
subplot(2,1,1); %分割窗口
plot(abs(yf)); %绘出原始信号的振幅图
axis([0, 12000 ,0,200]); %显示X轴、Y轴取值范围
title('滤波前的振幅图'); %加图形标题说明
yf1=fft(ys1); %对滤波后的声音进行FFT变换
subplot(2,1,2); %分割窗口
plot(abs(yf1)); %绘出滤波后信号的振幅图
axis([0, 12000 ,0,200]); %显示X轴、Y轴取值范围
title('滤波后的振幅图'); %加图形标题说明
通过滤波前后两信号的比较，通过上面的振幅图可以很清晰地观察到经过滤波后与滤波前存在明显差别，滤波后的图比较平稳，原始语音信号比滤波后的信号清晰，滤波后的声音显得比较低沉，因为原始信号经过低通滤波器后，低通滤波器已经把信号高频部分滤掉了，只剩下低频部分。
结论
经过一个星期的设计，本人基本按质按量完成了课程设计的各项工作。通过本次课程设计，使我加深了对数字信号处理这门课和MATLAB这个软件工具的认识。语音信号处理是以MATLAB为工具平台、数字信号处理为基础而形成的一门综合性学科，处理的目的是要得到一些语音参数以便高效的传输或存储，或者通过处理的某种运算以达到某种用途的要求。而MATLAB是一款多功能多特点的的软件，对于现代科学工程的计算和分析具有重要的意义。语音信号处理与滤波的设计对于所学的数字信号处理课程的理论得到了综合的应用，并从实践上初步实现对数字信号的处理。
在设计过程中遇到了很多从未遇到过的问题，小到一些标点符号，大到程序的设计，但我还是有耐心，经过多方面查资料和请教老师、同学，逐渐把一些不懂的问题解决，最终完成本次课程设计。在本次设计过程中，使我对数字信号处理的设计和MATLAB软件的认识更进一层，更重要的是使我学会了如何面对困难。在设计中，总结了以下几点经验，这也是自己在实践中一步一步摸索出来。
1、开始设计时，不要太过了想象完美，有思路就马上记录下来，在查阅资料也要注意记录重点。
2、遇到问题，不要着急，保住沉着的心态，通过查阅资料、请教老师或同学找出解决途径。
3、做设计，单凭自己的力量是不够的，要学会学习别人好的方面，多查阅资料，多与老师、同学交流，把别人的知识转变为自己的知识。
4、养成良好的工作态度，对以后工作具有具有重要意义。
在完成设计的过程中，由于所学知识的狭窄，本次设计不免有很多不理想的地方，这些都是有待提高的。
参考文献
(1)程佩青.数字信号处理教程(第三版).清华大学出版社.20##年2月
(2)赵红怡,张常年.数字信号处理及其MATLAB实现.化学工业出版社.20##年1月
(3)钱同惠.数字信号处理.机械工业出版社.20##年2月
(4)李海涛,邓樱.MATLAB6.1基础及应用技巧国防工业出版社20##年3月