数字信号处理期末综合实验报告.doc

数字信号处理综合实验报告实验题目：基于Matlab的语音信号去噪及仿真专业名称：学 号：姓 名：日 期：报告内容：一、 实验原理1、 去噪的原理1.1 采样定理在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中，最高频率fmax的2倍时，即：fs.max=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的510倍；采样定理又称奈奎斯特定理。 1924年奈奎斯特(Nyquist)就推导出在理想低通信道的最高大码元传输速率的公式: 理想低通信道的最高大码元传输速率=2W*log2 N (其中W是理想低通信道的带宽,N是电平强度)为什么把采样频率设为8kHz?在数字通信中，根据采样定理, 最小采样频率为语音信号最高频率的2倍频带为F的连续信号 f(t)可用一系列离散的采样值f(t1),f(t1t)，f(t12t)，.来表示,只要这些采样点的时间间隔t1/2F，便可根据各采样值完全恢复原来的信号f(t)。 这是时域采样定理的一种表述方式。 时域采样定理的另一种表述方式是：当时间信号函数f(t)的最高频率分量为fM时,f(t)的值可由一系列采样间隔小于或等于1/2fM的采样值来确定,即采样点的重复频率f2fM。图为模拟信号和采样样本的示意图。 时域采样定理是采样误差理论、随机变量采样理论和多变量采样理论的基础。对于时间上受限制的连续信号f(t)（即当tT时,f(t)=0,这里T=T2-T1是信号的持续时间），若其频谱为F（）,则可在频域上用一系列离散的采样值 （1-1）采样值来表示,只要这些采样点的频率间隔 （1-2）。1.2 采样频率采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz）来表示。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率只能用于周期性采样的采样器，对于非周期性采样的采样器没有规则限制。 采样频率的常用的表示符号是 fs。 通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据采样定理，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精确。对于采样率来说你可以想象它类似于一个照相机，44.1kHz意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样率越高，计算机摄取的图片越多，对于原始音频的还原也越加精确2、 数字滤波器设计的基本原理1.1 FIR数字滤波器的设计FIR：有限脉冲响应滤波器。有限说明其脉冲响应是有限的。与IIR相比，它具有线性相位、容易设计的优点。这也就说明，IIR滤波器具有相位不线性，不容易设计的缺点。而另一方面，IIR却拥有FIR所不具有的缺点，那就是设计同样参数的滤波器，FIR比IIR需要更多的参数。这也就说明，要增加DSP的计算量。DSP需要更多的计算时间，对DSP的实时性有影响。FIR滤波器的设计比较简单，就是要设计一个数字滤波器去逼近一个理想的低通滤波器。通常这个理想的低通滤波器在频域上是一个矩形窗。根据傅里叶变换我们可以知道，此函数在时域上是一个采样函数。通常此函数的表达式为：sa（n）sin（n）/n（1-3）但是这个采样序列是无限的，计算机是无法对它进行计算的。故我们需要对此采样函数进行截断处理。也就是加一个窗函数。就是传说中的加窗。也就是把这个时域采样序列去乘一个窗函数，就把这个无限的时域采样序列截成了有限个序列值。但是加窗后对此采样序列的频域也产生了影响：此时的频域便不在是一个理想的矩形窗，而是成了一个有过渡带，阻带有波动的低通滤波器。通常根据所加的窗函数的不同，对采样信号加窗后，在频域所得的低通滤波器的阻带衰减也不同。通常我们就是根据此阻带衰减去选择一个合适的窗函数。如矩形窗、汉宁窗、汉明窗、BLACKMAN窗、凯撒窗等。1.2 IIR数字滤波器的设计对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：（1）确定所需类型数字滤波器的技术指标。（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为=2/T tan(0.5) （1-4）(3)将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。（4）设计模拟低通滤波器。（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。二、 实验步骤 其大概流程框图可如下表示：（图2-1）语音信号采集效果显示、对比语音信号录入语音信号变换信号加噪语音信号滤波 图2-11. 滤波器的设计。 利用窗函数法设计FIR滤波器的步骤。如下：（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。（2）构造希望逼近的频率响应函数。（3）计算h(n).。（4）加窗得到设计结果。接下来，我们根据语音信号的特点给出有关滤波器的技术指标：低通滤波器的性能指标：fp=1000Hz，fc=1200Hz，As=50db ,Ap=1dB高通滤波器的性能指标：fp=3500Hz，fc=4000Hz，As=50dB，Ap=1dB在Matlab中，可以利用函数fir1设计FIR滤波器，利用Matlab中的函数freqz画出各步步器的频率响应。MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：B,A=cheby1(N,Rp,wpo,ftypr)B,A=cheby1(N,Rp,wpo,ftypr,s)利用模拟滤波器设计IIR数字低通滤波器的步骤：如下：（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。（3）按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：B,A=cheby1(N,Rp,wpo,ftypr)B,A=cheby1(N,Rp,wpo,ftypr,s)函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。2. 语音信号的录制。单击自己的电脑开始程序，选择所有程序，接着选择附件，再选择录音。自己录入语音信号，然后保存在MATLAB文件夹里面，命名为“chenghaijie1.wav”。3. 在MATLAB平台上读入语音信号、绘制频谱图并回放原始语音信号。利用MATLAB中的wavread命令来读入（采集）语音信号，将它赋值给某一向量。y,fs,bits=wavread( N1 N2);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。N1 N2表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值。 4. 利用matlab函数randn编程加入一段随机噪音信号，再利用设计的FIR和IIR滤波器去噪，分别绘制去噪后的频谱图、回放语音信号与原始信号的频谱图、原始语音信号比较，并且比较两种滤波器的优缺点和得出语音信号的频段。三、 实验结果下面我们将给出设计FIR数字滤波器的主要程序和图像。 FIR低通滤波器程序见附录1； FIR低通滤波器图像：（图22）图22 FIR低通滤波器 FIR高通滤波程序见附录2； FIR高通滤波图像：（图23）图23 FIR高通滤波器下面我们将给出IIR数字滤波器的主要程序。IIR低通滤波器程序见附录3；IIR低通滤波器图像：（图24）图24 IIR低通滤波器IIR滤波器高通程序见附录4；IIR滤波器高通图像：（图25）图25 IIR高通滤波器原始语音信号的回放及频谱分析程序：见附录5；原始语音信号及其频谱分析图像、运行程序可以听到原始信号的回放：（图26） 图26 原始语音信号及其频谱分析图像利用randn函数把一段随机噪音信号加入原始语音信号的信号处理过程程序：见附录6；加噪后语音信号的时域波形、频谱图：（图27） （图27）加噪后语音信号的时域波形、频谱图利用FIR低通滤波器法去噪程序：见附录7；FIR滤波器去噪前与去噪后语音信号的时域波形、频谱图：（图28） （图28）FIR滤波器去噪前与去噪后语音信号的时域波形、频谱图利用IIR低通滤波器法去噪程序：见附录8；IIR滤波器去噪前与去噪后语音信号的时域波形、频谱图：（图29） （图29）IIR滤波器去噪前与去噪后语音信号的时域波形、频谱图四、 实验分析1. 由（图28）FIR滤波器去噪前与去噪后语音信号的时域波形、频谱图和（图29）IIR滤波器去噪前与去噪后语音信号的时域波形、频谱图可看出原始语音信号和加噪语音信号时域波形和频谱图的区别。加噪后的语音信号的时域波形比原始语音信号要模糊得多，频谱图则是在频率5000Hz以后出现了明显的变化。再通过滤波前的信号波形和频谱图的对比，可以明显看出滤波后的波形开始变得清晰了，有点接近原始信号的波形图了。滤波后信号的频谱图也在5000Hz以后开始逐渐接近原始语音信号的频谱图。再从对语音信号的回放，人耳可以明显辨别出加噪后的语音信号比较浑浊，还有很明显嘎吱嘎吱的杂音在里面。滤波后，语音信号较加噪后的信号有了明显的改善，基本可以听清楚了，而且杂音也没有那么强烈，但是声音依然没有原始语音信号那么清晰脆耳。2. 结合去噪后的频谱图对比可知FIR滤波器和IIR滤波器的优缺点IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。 （1）IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。FIR滤波器则要求较低。（2）IIR滤波器运算误差大，有可能出现极限环振荡，FIR相比之下运算误差较小，不会出现极限环振荡。（3）IIR幅频特性精度很高，不是线性相位的，可以应用于对相位信息不敏感的音频信号上； （4）与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。 （5）IIR单位响应为无限脉冲序列FIR单位响应为有限的 （6）FIR幅频特性精度较之于iir低，但是线性相位，就是不同频率分量的信号经过FIR滤波器后他们的时间差不变。这是很好的性质。 （7）IIR滤波器有噪声反馈，而且噪声较大，FIR滤波器噪声较小。FIR幅频特性精度较之于IIR低，但是线性相位，就是不同频率分量的信号经过FIR滤波器后他们的时间差不变，这是很好的性质。3. 通过对语音信号做频谱分析和通过比较加噪前后，语音的频谱和语音回放，能明显的感觉到加入噪声后回放的声音与原始的语音信号有很大的不同,前者随较尖锐的干扰啸叫声。从含噪语音信号的频谱图中可以看出含噪声的语音信号频谱,在整个频域范围内分是布均匀。其实，这正是干扰所造成的。通过滤波前后的对比，低通滤波后效果最好，高通滤波后的效果最差。由此可见，语音信号主要分布在低频段，而噪声主要分布在高频段。五、 实验总结基于Matlab的语音信号去噪及仿真课题的特色在于它将语音信号看作一个向量，于是就把语音数字化了。那么，就可以完全利用数字信号处理的知识来解决语音及加噪处理问题。我们可以像给一般信号做频谱分析一样，来对语音信号做频谱分析，也可以较容易的用数字滤波器来对语音进行滤波处理。通过比较加噪前后，语音的频谱和语音回放，能明显的感觉到加入噪声后回放的声音与原始的语音信号有很大的不同,前者随较尖锐的干扰啸叫声。从含噪语音信号的频谱图中可以看出含噪声的语音信号频谱,在整个频域范围内分是布均匀。其实，这正是干扰所造成的。通过两种滤波器滤波前后的对比，可知FIR滤波器和IIR滤波器的优缺点，低通滤波后效果最好，高通滤波后的效果最差。由此可见，语音信号主要分布在低频段，而噪声主要分布在高频段。六、 参考文献1 Sanjit K.Mitra，数字信号处理实验指导书，电子工业出版社，2005 年 1月 2 胡航，语音信号处理，哈尔滨工业大学出版社，2000 年 5 月 3姚天任.数字语音处理M.武汉：华中科技大学出版社，2005 年 3 月七、 附录附录 1Ft=8000;%设置抽样频率为8000HzFp=1000;%设置通带边界为1000HzFs=1200;%设置阻带边界为1200Hzwp=2*Fp/Ft;%计算归一化通带边界角频率ws=2*Fs/Ft;%计算归一化阻带边界角频率rp=1;%设置通带纹波为1dBrs=50;%设置阻带衰减为50dBp=1-10.(-rp/20);%通带阻带波纹s=10.(-rs/20);fpts=wp ws;mag=1 0;dev=p s;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率b21=fir1(n21,wn21,Kaiser(n21+1,beta);%由fir1设计滤波器h,w=freqz(b21,1);%得到频率响应plot(w/pi,abs(h);%绘制FIR低通滤波器title(FIR低通滤波器);%标题附录2Ft=8001;%设置抽样频率为8001HzFp=4000;%设置通带边界为4000HzFs=3500;%设置阻带边界为3500Hzwp=2*Fp/Ft;%计算归一化通带边界角频率ws=2*Fs/Ft;%计算归一化阻带边界角频率rp=1;%设置通带纹波为1dBrs=50;%设置阻带衰减为50dBp=1-10.(-rp/20);%通带阻带波纹s=10.(-rs/20);fpts=ws wp;mag=0 1;dev=p s;n23,wn23,beta,ftype=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率b23=fir1(n23,wn23,high,Kaiser(n23+1,beta);%由fir1设计滤波器h,w=freqz(b23,1);%得到频率响应plot(w*12000*0.5/pi,abs(h);%绘制FIR高通滤波器title(FIR高通滤波器);%标题axis(3000 6000 0 1.2);%确定横纵坐标的范围附录3Ft=8000;%设置抽样频率为8000HzFp=1000;%设置通带边界为1000HzFs=1200;%设置阻带边界为1200Hzwp=2*pi*Fp/Ft;%计算归一化通带边界角频率ws=2*pi*Fs/Ft;%计算归一化阻带边界角频率fp=2*Ft*tan(wp/2);%计算通带边界频率fs=2*Fs*tan(wp/2);%计算阻带边界频率n11,wn11=buttord(wp,ws,1,50,s);%求低通滤波器的阶数和截止频率b11,a11=butter(n11,wn11,s);%利用butter求S域的频率响应的参数num11,den11=bilinear(b11,a11,0.5);%双线性变换实现S域到Z域的变换h,w=freqz(num11,den11);%根据参数求出频率响应plot(w*8000*0.5/pi,abs(h);%绘制IIR低通滤波器title(IIR低通滤波器);%标题legend(用butter设计);%添加图例的标注grid;%使绘制的图有网格附录4Ft=8001;%设置抽样频率为8001HzFp=4000;%设置通带边界为4000HzFs=3500;%设置阻带边界为3500Hzwp1=tan(pi*Fp/Ft);%通带到低通滤波器参数转换ws1=tan(pi*Fs/Ft);%阻带到低通滤波器参数转换wp=1;%计算归一化通带边界角频率ws=wp1*wp/ws1;%计算归一化阻带边界角频率n13,wn13=cheb1ord(wp,ws,1,50,s);%求模拟的低通滤波器阶数和截止频b13,a13=cheby1(n13,1,wn13,s);%利用cheby1求S域的频率响应的参数num,den=lp2hp(b13,a13,wn13); %将S域低通参数转为高通的num13,den13=bilinear(num,den,0.5); %利用双线性变换实现S域到Z域转h,w=freqz(num13,den13);%得到频率响应plot(w*21000*0.5/pi,abs(h);%绘制IIR高通滤波器title(IIR高通滤波器);%标题legend(用cheby1设计);%添加图例的标注附录5x,fs,bits=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令来读入（采集）语音信号，将它赋值给某一向量sound(x,fs,bits);%播放语音信号X=fft(x,100000);%对信号做50000点FFT变换magX=abs(X);%求幅度angX=angle(X);%求相位plot(x);title(原始信号波形);%绘制原始信号的波形plot(X); title(原始语音信号采样后的频谱图)%绘制原始语音信号采样后的频谱图附录6y,fs,bits=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令来读入（采集）语音信号，将它赋值给某一向量sound(y,fs)%播放语音信号n=length(y)%计算向量y的长度y_p=fft(y,n);%对向量y进行fft变换f=fs*(0:n/2-1)/n;%对应点的频率figure(1)%建立第一幅图subplot(2,1,1);%分配绘制的位置为第一块位置plot(y);%绘制原始语音信号采样后时域波形title(原始语音信号采样后时域波形);%标题xlabel(时间轴)%标注X轴ylabel(幅值 A)%标注Y轴subplot(2,1,2);%分配绘制的位置为第二块位置plot(f,abs(y_p(1:n/2);%绘制原始语音信号采样后的频谱图title(原始语音信号采样后的频谱图);xlabel(频率Hz);%标注X轴ylabel(频率幅值);%标注Y轴L=length(y)%计算向量y的长度noise=0.1*randn(L,2);%设置加入的随机噪声函数y_z=y+noise;%将原始信号与随机噪声函数相叠加sound(y_z,fs)%试听加入随机噪声函数的信号n=length(y);%计算向量y的长度y_zp=fft(y_z,L);%对加入正弦噪声向量y进行fft变换f=fs*(0:L/2-1)/L;%计算对应点的频率figure(2)%建立第二幅图subplot(2,1,1);%分配绘制的位置为第一块位置plot(y_z);%绘制加噪语音信号时域波形title(加噪语音信号时域波形);%标题xlabel(时间轴)%标注X轴ylabel(幅值 A)%标注Y轴subplot(2,1,2);%分配绘制的位置为第二块位置plot(f,abs(y_zp(1:L/2);%绘制加噪语音信号频谱图title(加噪语音信号频谱图);%标题xlabel(频率Hz);%标注X轴ylabel(频率幅值);%标注Y轴附录7y,fs,bits=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令来读入（采集）语音信号，将它赋值给某一向量sound(y,fs)%播放语音信号n=length(y)%计算向量y的长度y_p=fft(y,n);%对向量y进行fft变换f=fs*(0:n/2-1)/n;%对应点的频率figure(1)%建立第一幅图subplot(2,1,1);%分配绘制的位置为第一块位置plot(y);%绘制原始语音信号采样后时域波形title(原始语音信号采样后时域波形);%标题xlabel(时间轴)%标注X轴ylabel(幅值 A)%标注Y轴subplot(2,1,2);%分配绘制的位置为第二块位置plot(f,abs(y_p(1:n/2);%绘制原始语音信号采样后的频谱图title(原始语音信号采样后的频谱图);xlabel(频率Hz);%标注X轴ylabel(频率幅值);%标注Y轴L=length(y)%计算向量y的长度noise=0.1*randn(L,2);%设置加入的随机噪声函数y_z=y+noise;%将原始信号与随机噪声函数相叠加sound(y_z,fs)%试听加入随机噪声函数的信号n=length(y);%计算向量y的长度y_zp=fft(y_z,L);%对加入正弦噪声向量y进行fft变换f=fs*(0:L/2-1)/L;%计算对应点的频率figure(2)%建立第二幅图subplot(2,1,1);%分配绘制的位置为第一块位置plot(y_z);%绘制加噪语音信号时域波形title(加噪语音信号时域波形);%标题xlabel(时间轴)%标注X轴ylabel(幅值 A)%标注Y轴subplot(2,1,2);%分配绘制的位置为第二块位置plot(f,abs(y_zp(1:L/2);%绘制加噪语音信号频谱图title(加噪语音信号频谱图);%标题xlabel(频率Hz);%标注X轴ylabel(频率幅值);%标注Y轴Ft=5000;%设置抽样频率为5000HzFp=1000;%设置通带边界为1000HzFs=1200;%设置阻带边界为1200Hzwp=2*Fp/Ft;%计算归一化通带边界角频率ws=2*Fs/Ft;%计算归一化阻带边界角频率rp=1;%设置通带纹波为1dBrs=50;%设置阻带衰减为50dBp=1-10.(-rp/20);%通带阻带波纹s=10.(-rs/20);fpts=wp ws;mag=1 0;dev=p s;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev);%kaiserord求阶数截止频率b21=fir1(n21,wn21,Kaiser(n21+1,beta);%由fir1设计滤波器h,w=freqz(b21,1);%得到频率响应plot(w/pi,abs(h);%绘制FIR低通滤波器title(FIR低通滤波器);%标题x=fftfilt(b21,y_z);%重叠相加法实现线性卷积的计算X=fft(x,n);%计算x的fftfigure(3);%建立第三幅图subplot(2,2,1);plot(f,abs(y_zp(1:n/2);%绘制滤波前信号的频谱title(滤波前信号的频谱);%标题subplot(2,2,2);plot(f,abs(X(1:n/2);%绘制滤波后信号的频谱title(滤波后信号的频谱);%标题subplot(2,2,3);plot