回声信号检测系统设计28终期报告(共17页)

1、信号与系统A（2）课程(kchng)项目实施(shsh)报告 题目(tm)( B )： 回声信号检测系统设计 组 号： B28 任课教师： 组 长：_ _成 员：_成 员： _成 员： _成 员：_联系方式： 二零 15年 5 月 22 日目录(ml)：一、课程(kchng)项目实施方案二、项目的设计(shj)过程单次回声的分析与消除1、回声的产生及信号的时域、频域分析2、自相关函数分析回声延时量3、用迭代法求回声衰减系数4、回声的消除多次回声的分析，与消除1、回声的产生和回声延时量的分析2、用归一化法估算回声衰减系数3、回声的消除三、项目的结果分析四、分工合作情况参

2、考文献附录（matlab代码）课程(kchng)项目实施方案题目(tm)要求：B组 回声信号检测系统(xtng)设计设计一个回声信号检测系统，能够检测出原信号中包含的回声信号，并分析回声信号的参数。要求：（1）输入混有回声的声音（或信号源）；（2）设计检测系统，分析回声信号产生的原因、分析回声信号的幅度、分析回声信号的延迟量； (3) 能够现场演示。设计思想：我们组的设计思想是利用matlab软件平台，用程序人为制造回声，并且在仅知道回声信号的条件下，分析延时量和衰减系数，从而消除回声，并且测出声源与障碍物的距离。项目的设计过程单次回声的分析与消除1、回声的产生及信号的时域、频域分析对于已经采

3、集好的一段语音信号，利用matlab的”audioread”命令可以对该信号进行读取和采样，本项目中采用的文件为wav文件，采样频率为44100.x,fs=audioread(E:ceshi.wav);利用plot和fft函数可以对已采集的信号分别绘制时域波形图和经快速傅里叶变换得到的频谱图；abs和angle函数可以绘制信号频域的幅值与相位图。 利用zeros函数创建空矩阵可以对原始(yunsh)信号进行延时，并产生回声。本例子中，假设(jish)延时量为8000个采样点，衰减系数为0.4。同样对回声(hushng)信号y进行时域、频域分析。2、自相关函数分析(fnx)回声延时量自相关函数(

4、hnsh)的数学定义为：同一时间函数在瞬时t和t+a的两个(lin )值相乘积的平 HYPERLINK /view/1052684.htm t _blank 均值作为延迟时间t的函数，它是信号与延迟后信号之间相似性的度量。简单地说，自相关函数是用来表征一个信号在不同时刻t1,t2状态之间的相关程度，其在某一时间点的值越大，相关程度越高。在matlab中，对于一个混有回声的语音信号，可以利用xcorr自相关函数可以求出回声信号与原始信号的采样点之差，进而求得延时量。 从自相关函数图像中可以看到，在采样点的中心(zhngxn)有一最大的幅值，其含义就是当延时量为零时，信号相关程度最高；而在这一峰值

5、的两侧有两个次大的幅值，其含义就是在该处信号相关程度次高；而在别的点处其相关程度相对于该2个峰值都可以忽略，因此算出这2个峰值所对应的采样点之差，就是回声信号中，原始信号与回声信号对应的采样点之差。由此，可以计算(j sun)出该信号中采样点之差为8000个单位，这与我们所假设的值是吻合的。最后，将这个采样点之差除以采样频率，即，便可得到(d do)实际的延时量。3、用迭代法求回声衰减系数已知回声的延迟量后，还需要知道回声的衰减系数才可以确定回声信号。对于一个一次回声信号，我们可以采用迭代的方法来求衰减系数:设xn为原始信号，yn为原始信号与回声信号的叠加则有:y(1)=x(1)Y(1+N)=

6、X(1+N)+aX(1)Y(1+2N)=x(1+2N)+ax(1+N)y(1+(k-1)N)=x(1+(k-1)N)+ax(1+(k-2)N)Y(1+kN)=ax(1+(k-1)N)上述迭代(di di)的公式可以由下面线段图来形象(xngxing)的说明(shumng)：y序列的第一个点是没有回声的，在第（1+N）、（1+2N）、（1+（k-1）N）的点都是由原始信号与衰减后回声信号叠加而成的，最后一个（1+kN）个点是仅含衰减后的回声的。由图像可以看出:设y(n)的长度为L，对于k,则满足1+kNL（且k是整数）这个k是唯一的并且可求的。因此由最后一个式子逐级带回上面的式子，将含x(n)的

7、项全迭代成含y(n)和衰减系数的项，最后可以得到一个关于的一元k次方程。1/a*y(1+kN)=y(1+(k-1)N)-ay(1+(k-2)N)+a2*y(1+(k-3)N)-a3*y(1+(k-4)N)+.ak*y(1)接着用matlab中的roots函数可以对进行求解。但是对于求解出的k个根，排除虚数和不在01范围内的数，可能仍有一些根是满足条件的。这时候我们可以将上述的迭代公式从y(2)开始计算，即：y(2)=x(2)Y(2+N)=X(2+N)+aX(2)Y(2+2N)=x(2+2N)+ax(2+N)y(2+(k-1)N)=x(2+(k-1)N)+ax(2+(k-2)N)Y(2+kN)=

8、ax(2+(k-1)N)对于同一个信号衰减系数必定相同，那么在两个方程中，相同的根便是所求的。我们组为了实现(shxin)程序的的普遍性，对程序进行了优化，自动(zdng)找出两组解中相同的解。因此(ync)用迭代法可以精准的算出回声的衰减系数。4、回声的消除由上面求解(qi ji)得到的延时量N_delay=8000，衰减系数=0.4，可以(ky)写出原始信号与含回声信号的关系式y(n)=x(n)+0.4x(n-8000),其系统函数为H(z)= QUOTE 在matlab中可以用filter函数对y信号(xnho)进行滤波：z即为滤波完后的信号，sound(z,fs)发现单次回声已经全部消

9、除。多次回声的分析，与消除1、回声的产生和回声延时量的分析多次回声的产生原理与一次回声类似，利用zeros函数创建空矩阵可以对原始信号进行延时。本例中延时量分别为8000个和12000个采样点，衰减系数分别为0.4和0.2.同样的，用xcorr自相关函数计算延迟量。 其中(qzhng)N_delay1是第一次回声延迟的采样点，N_delay2是第二次回声延迟的采样点。2、用归一化法估算( sun)回声衰减系数加入(jir)两次回声后，无法采用之前(zhqin)单次回声的方法计算衰减系数，因为无法确定唯一的k来列出迭代的方程式。我们采用自相关函数的归一化法估算回声衰减系数。具体过程如下：设Ry(

10、n)为y(n)序列的自相关函数，Rx(n)为x(n)序列的自相关函数。由自相关函数的定义可以化简得到以下方程：Ry(n)=(1+ QUOTE + QUOTE )Rx(n) + QUOTE Rx(n- QUOTE ) +Rx(n+ QUOTE ) QUOTE Rx(n- QUOTE ) +Rx(n+ QUOTE )+ QUOTE Rx(n+ QUOTE - QUOTE ) +Rx(n+ QUOTE )将n=0, QUOTE ,N2分别代入R-1(0)=(1+ QUOTE + QUOTE )Rx(0)+ QUOTE Rx(- QUOTE ) +Rx( QUOTE ) QUOTE Rx(- QUOT

11、E ) +Rx( QUOTE )+ QUOTE Rx( QUOTE - QUOTE ) +Rx( QUOTE ) （2）考虑到Rx(- QUOTE ) =Rx( QUOTE )约等于0，Rx(- QUOTE ) =Rx( QUOTE )约等于0，Rx( QUOTE - QUOTE ) =Rx( QUOTE )约等于0，（2）式化为Ry (0)=(1+ QUOTE + QUOTE )Rx(0) （3）同理有Ry( QUOTE )约等于 QUOTE Rx(0) ，Ry( QUOTE )约等于 QUOTE Rx(0)因此有 QUOTE , QUOTE 还可推广到n次回声由此可以(ky)估算出第K次衰

12、减系数。 利用(lyng)matlab的fsolve解非线性方程组的功能，便可估算( sun)出2个回声的衰减量。与理论的0.4和0.2有5%左右的误差。3、回声的消除同样利用filter函数。H(z)= QUOTE 执行sound(z,fs)后，回声明显不见了，可以达到效果。项目的结果分析利用回声测距上述(shngsh)两例单次回声和多次回声的消除效果(xiogu)是良好的，消除后的声音完全听不出回声。此外通过(tnggu)对回声的分析可以应用与测距。其原理很简单：通过自相关函数我们已求出回声的延迟采样点个数，N_delay。其真实时间其中Fs是采样率，t的单位为秒。则声源与障碍物的距离。以

13、报告中的例子举例：单次回声中N_delay=8000，Fs=44100,则可算出t=0.18s，那么声源与障碍物的距离s=340*0.18/2=30.6m四、分工合作情况查找资料：盛亦池 胡斌 刘杨雨编写代码：盛亦池 应祺超代码调试：胡斌 应祺超 马嘉敏终期报告：马嘉敏 刘杨雨参考文献：1、自相关函数 百度百科/link?url=OWeCZQXASKxx4INScryWy5uXYERPm00lFVNbzuRY-H-s8LgSyr63aq4-Pej7Vi9C3YjS4ah6jyhJjEuOgtykwa2、matlab课程设计-基于MATLAB的回波信号的产生与消除 道客巴巴 HYPERLINK

14、/p-7062907417869.html /p-7062907417869.html3、信号与系统实验指导/p-4147988583065.html4、基于MATLAB的数字信号处理与实践开发，王佳梅，西安电子科技大学出版社附录一次回声matlab代码：%制造回声x,fs=audioread(E:ceshi.wav);N=length(x);x1=x(1:N);x2=x(1:N);x1=x1,zeros(1,8000);x2=zeros(1,8000),0.4*x2;y=x1+x2;%原始信号(xnho)时域、频域、幅值、相位图n1=0:22050;figure(1)subplot(2,2,

15、1);plot(x);title(原始(yunsh)信号时域图);X=fft(x);x1=x(1:4:length(x);X1=x1(:,1);df=fs/length(X);Fx=df*(0:length(X)-1);subplot(2,2,2);plot(Fx(1:44100/4),X(1:44100/4);title(原始(yunsh)信号频谱图);subplot(2,2,3);plot(abs(X);title(原始信号幅值);subplot(2,2,4);plot(angle(X);title(原始信号相位);%叠加信号的时域、频域、幅值、相位图figure(2);subplot(2

16、,2,1);plot(y);title(带有回声回声时域图);Y=fft(y);y1=y(1:4:length(y);Y1=y1(:,1);df2=fs/length(Y);Fy=df2*(0:length(Y)-1);subplot(2,2,2);plot(Fy(1:44100/4),Y(1:44100/4);title(带有回声信号频谱图);subplot(2,2,3);plot(abs(Y);title(回声幅值);subplot(2,2,4);plot(angle(Y);title(回声(hushng)相位);%自相关(xinggun)函数 算延时量cy=xcorr(y);figure

17、(3)plot(cy);title(多次回声的自相关函数(hnsh)图像);h1_y,h1=max(cy);cy2=cy;cy2(1,h1-1000:h1+1000)=0;h2_y,h2=max(cy2);N_delay=abs(h1-h2);kmax=floor(length(y)+N_delay-1)/N_delay);for k=1:kmax t1(k)=(-1)k*y(1,N_delay*(k-1)+1);endjie1=roots(t1);for i=1:length(jie1) if isreal(jie1(i) jie1(i)=jie1(i); else jie1(i)=0; e

18、ndendfor k=1:kmax t2(k)=(-1)k*y(1,N_delay*(k-1)+2);endjie2=roots(t2);for i=1:length(jie2) if isreal(jie2(i) jie2(i)=jie2(i); else jie2(i)=0; endendjie1(find(jie1(:,1)=0)= ;jie2(find(jie2(:,1)=0)= ;jie1=jie1;jie2=jie2;jie=(intersect(int16(jie1*1000),int16(jie2*1000);format shortfinaljie=double(jie)/1

19、000;%消除(xioch)回声a=zeros(1,N);a(1)=1;a(N_delay+1)=finaljie;z=filter(1,a,y);多次回声(hushng)matlab代码：x,fs=audioread(E:ceshi.wav);%制造(zhzo)回声N=length(x);x1=x(1:N);x2=x(1:N);x3=x(1:N);x1=x1,zeros(1,12000);x2=zeros(1,8000),0.4*x2,zeros(1,4000);x3=zeros(1,12000),0.2*x3;y=x1+x2+x3;%原始信号时域、频域、幅值、相位图n1=0:22050;f

20、igure(1)subplot(2,2,1);plot(x);title(原始信号时域图);X=fft(x);x1=x(1:4:length(x);X1=x1(:,1);df=fs/length(X);Fx=df*(0:length(X)-1);subplot(2,2,2);plot(Fx(1:44100/4),X(1:44100/4);title(原始信号频谱图);subplot(2,2,3);plot(abs(X);title(原始信号幅值);subplot(2,2,4);plot(angle(X);title(原始相位);%叠加信号的时域、频域、幅值、相位图figure(2);subplot(2,2,1);plot(y);title(带有回声信号时域图);Y=fft(y);y1=y(1:4:length(y);Y1=y1(:,1);df2=fs/length(Y);Fy=df2*(0:length(Y)-1);subplot(2,2,2);plot(Fy(1:44100/4),Y(1:44100