MATLAB实现数字信号处理

《数字信号处理》课程设计实例：声音信号的处理一．摘要：这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

通过对声音的采样，将声音采样后的频谱与滤波。

MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。

经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。

MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。

在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。

二．课程设计目的：综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

三．设计容：容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音容，分析两段容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。

四．设计原理：4.1．语音信号的采集熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有：a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。

Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据；b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y；c：wavwrite((y,fs,wavfile);创建音频文件；d：y=wavread(file);读取音频文件；关于声音的函数还有sound()；soundsc()；等。

4.2滤波器：4．21.IIR滤波器原理冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。

为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。

双线性法设计IIR 数字滤波器的步骤： 1）将数字滤波器的频率指标{ kΩ}由Wk ＝（2/T ）*tan(kw ),转换为模拟滤波器的频率指标{kΩ}.2)由模拟滤波器的指标设计H(s).3)由H(s)转换为H(z)1121()()1z H z H s s T z ---==+4.22.FIR 滤波器原理FIR 滤波器与IIR 滤波器特点不同，IIR 滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。

而有限长单位冲激响应（FIR ）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。

由于FIR 系统的冲激响应就是其系统函数各次项的系数，所以设计FIR 滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H （z ）的系数，冲激响应长度N 就是系统函数H （z ）的阶数。

只要N 足够长，截取的方法合理，总能满足频域的要求。

这种时域设计、频域检验的方法一般要反复几个回合，不像IIR DF 设计靠解析公式一次计算成功。

给出的理想滤波器频率响应是，它是w 的周期函数，周期，因此可展开成傅氏级数由傅立叶反变换导出，即1()()2jw jwnd d h n He e dw π=⎰，再将()d h n 与窗函数()w n 相乘就可以得到()h n ，。

、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。

但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。

窗函数的傅式变换W(ej ω)的主瓣决定了H(ej ω)过渡带宽。

W(ej ω)的旁瓣大小和多少决定了H(ej ω)在通带和阻带围波动幅度，常用的几种窗函数有：矩形窗 w(n)=RN(n)；Hanning 窗；Hamming窗；Blackmen窗；Kaiser窗。

式中Io(x)为零阶贝塞尔函数。

五．设计步骤：5.1录制女音：利用MATLAB中的函数录制声音。

function nvyin()fs=11025; %采样频率str=[''];disp(str);disp(' 开始录音');str=[''];disp(str);y=wavrecord(3*fs,fs,'double'); %录制声音3秒str=[''];disp(str);disp(' 录音结束');str=[''];disp(str);str=[''];disp(str);disp(' 播放录音');str=[''];disp(str);wavplay(y,fs); %播放录音str=[''];disp(str);disp(' 播放录音结束');str=[''];disp(str);wavwrite(y,fs,'原女音'); %声音的存储5.2采样语音信号并画出时域波形和频谱图读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。

程序如下：[x,fs,bits]=wavread('女音.wav'); %读取声音N=length(x); %计数读取信号的点数t=(1:N)/fs; %信号的时域采样点f0=fs/N; %采样间隔n=1:N/2; %取信号的一半figure(1);subplot(2,2,1); %把画图区域划分为2行2列，指定第一个图plot(t, x); %画出声音采样后的时域波形title('原女音信号的时域波形'); %给图形加注标签说明xlabel('时间/t');ylabel('振幅/A');grid ; %添加网格y=fft(x); %对信号做N点FFT变换k=(n-1)*f0; %频域采样点subplot(2,2,3); %把画图区域划分为2行2列，指定第三个图plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图title('FFT变换后声音的频谱特性'); %给图形加注标签说明xlabel('频率/Hz');ylabel('幅值/A');grid; %添加网格subplot(2,2,4); %把画图区域划分为2行2列，指定第四个图if y~=0 %判断指数是否为0plot(k,20*log10(abs(y(n)))); %画信号频谱的分贝图endxlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid ; %添加网格%实际发出声音落后录制动作半拍的现象的解决siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]); %截取语音信号N=length(x1); %计数读取信号的点数t=(1:N)/fs; %信号的时域采样点f0=fs/N; %采样间隔n=1:N/2; %取信号的一半figure(2);subplot(2,2,1); %把画图区域划分为2行2列，指定第一个图plot(t,x1); %画出声音采样后的时域波形title('截取后女音信号的时域波形'); %给图形加注标签说明xlabel('时间/t');ylabel('振幅/A');grid ; %添加网格y1=fft(x1); %对信号做N点FFT变换subplot(2,2,3); %把画图区域划分为2行2列，指定第三个图k=(n-1)*f0; %频域采样点plot(k,abs(y(n))); %绘制原始语音信号的幅频响应图title('FFT变换后声音的频谱特性'); %给图形加注标签说明xlabel('频率/Hz');ylabel('幅值/A');grid ; %添加网格subplot(2,2,4); %把画图区域划分为1行2列，指定第二个图 if y1~=0 %判断指数是否为0 plot(k,20*log10(abs(y1(n)))); %画信号频谱的分贝图 end xlabel('Hz'); ylabel('振幅/分贝');title('FFT 变换后声音的频谱特性');grid ; %添加网格0123-1-0.500.51原女音信号的时域波形时间/t振幅/A200040006000100200300FFT 变换后声音的频谱特性频率/Hz幅值/A0200040006000-5050Hz振幅/分贝FFT 变换后声音的频谱特性0123-1-0.500.51截取后女音信号的时域波形时间/t振幅/A200040006000100200300FFT 变换后声音的频谱特性频率/Hz幅值/A0200040006000-5050Hz振幅/分贝FFT 变换后声音的频谱特性结果分析：由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。