matlab课程设计-基于MATLAB的回波信号的产生与消除

数字信号处理课程设计题目：基于MATLAB的回波信号的产生与消除课程：MATLAB课程设计姓名：学号：摘要在这个课程设计中，利用matlab采集一段语音，在这段语音的基础上，加入一定延时和衰减的回音，最后消去回音并且测出延时时间来计算障碍物距离正文①设计目的与要求采集语音：采集一段语音，绘制其时域波形，对此音频信号用FFT作谱分析。加入回声：对采集的语音进行处理，加入一段回声，并绘制其时域波形，对其进行FFT频谱分析，绘制频谱图。从带有回声的声音信号中恢复原信号：设计合适的滤波器，对带有回声的声音信号进行滤波，恢复原信号。绘制所设计滤波器的幅频和相频特性，及滤波后的信号的时域波形和频谱图。从带有回声的声音信号中估计反射物的距离：采用相关分析法从带有回声的声音信号中估计反射物的距离。②具体内容及原理（1）语音采集利用matlab采集一段语音并保存，代码如下fs=8000;x=wavrecord(3*fs,fs,'double');wavplay(x,fs);wavwrite(x,'原始信号');//存储音频：原始信号（2）原始信号的时域波形，FFT频谱分析代码如下subplot(3,1,1);plot(x);gridon;xlabel('时间');ylabel('幅值');title('原始信号时域波形');subplot(3,1,2);wx=fft(x);f=(0:3*fs-1)*fs/(3*fs);plot(f,abs(wx));gridon;xlabel('频率');ylabel('幅值');title('幅频特性');subplot(3,1,3);plot(f,angle(wx));gridon;xlabel('频率');ylabel('相位');title('相频特性');图如下：（3）加入回声在已有声音信号x的基础上产生带回声的声音信号，可以表达为在原信号的基础上叠加其延时衰减的分量。假设只有一个回声的情况下，可简化其模型为：y(n)=x(n)+ax(n-N)a为反射系数；N为延迟时间。在这里，取a=0.5，N=2400（即0.3秒的延时）下面则加入回声且保存代码如下：N=2400;y=[x;zeros(N,1)]+0.5*[zeros(N,1);x];wavwrite(y,'加回声后的信号');加回声后信号附件(双击打开):加回声后的信号.wav(4)加回声后信号的时域波形，FFT频谱分析代码如下：subplot(3,1,1);plot(y);gridon;xlabel('时间');ylabel('幅值');title('加回声后信号时域波形');subplot(3,1,2);wy=fft(y);f=(0:3*fs+N-1)*fs/(3*fs+N);plot(f,abs(wy));gridon;xlabel('频率');ylabel('幅值');title('幅频特性');subplot(3,1,3);plot(f,angle(wy));gridon;xlabel('频率');ylabel('相位');title('相频特性');图如下：(5)从带有回声的声音信号中恢复原信号且估计反射物的距离这里把信号的恢复和反射物距离的估计放到一起是基于这么一种考虑，说明如下：在回声产生的过程中，用到了：y(n)=x(n)+ax(n-N)，用的a=0.5,N=2400。然而现在要从加回声后的信号中恢复原信号，应该是在这么一种前提下，即“只有y(n)已知，其他都是未知的”。就是说，要假设我们并不知道原信号，且a与N都是未知的，这就给信号的恢复带来了困难，如果直接用y(n)=x(n)+0.5*x(n-2400)是不合理的。这个时候就要用到对反射物距离的估计的过程，在这个过程中利用相关分析法可以估算出N的值，利用N来算反射物的距离，求得N，则可以进一步求得a，具体方法和原理如下：如何求N利用自相关函数xcorr来估计N，对于信号x(n)，其长度为N，其求得的自相关函数为r(m)=1()*()Nnxnxnm，其中m的范围为-(N-1)到N-1，而且显然是左右对称的。下面是一个简单例子运行代码xcorr([123])结果如下ans=3.00008.000014.00008.00003.0000自相关函数是对函数本身在两个时刻t1，t2的相关程度的一种衡量标准，对于加回声后信号y(n)=x(n)+ax(n-N),y(n)是由x(n)与它的一个衰减延时ax(n-N)叠加而成，因为相关函数是函数本身相关程度的一种衡量，可以看到，y(n)的自相关函数将出现几个极值点，自相关函数为r(m)=1()*()Nnxnxnm，对于y(n)极值点应该出现在m=0，m=+-N，这时候相关程度相对较大。所以只要求出两个极值点之间的距离就能得到N如何求a知道N后该怎么求得a值，下面是一种想法：y(n)是x(n)前补零与后补零的叠加，于是有y(1)=x(1)y(1+N)=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)设y(n)的长度为L，对于k，则满足1+kNL,1+(k+1)NL现在逐级代入，可得1/a*y(1+kN)=x(1+(k-1)N)=y(1+(k-1)N)-ax(1+(k-2)N)=y(1+(k-1)N)-ay(1+(k-2)N)+a^2*x(1+(k-3)N)=...=y(1+(k-1)N)-ay(1+(k-2)N)+a^2*y(1+(k-3)N)-a^3*y(1+(k-4)N)+...a^k*y(1)即1/a*y(1+kN)=y(1+(k-1)N)-ay(1+(k-2)N)+a^2*y(1+(k-3)N)-a^3*y(1+(k-4)N)+...a^k*y(1)以此为方程来求解a为方便理解，下面是一个简单例子取x(n)=[12345]a=0.5N=2则有x(n)=12345x(n-2)=0.511.522.5y(n)=123.556.522.5有y(1)=x(1)y(3)=x(3)+ax(1)y(5)=x(3)+ax(5)y(7)=ax(5)得到1/a*y(7)=y(5)-ay(3)+a^2*y(1)即a^3*y(1)-a^2*y(3)+a*y(5)-y(7)=0把y(1)，y(3)，y(5)，y(7)代入，用matlab求解a，代码如下roots([1-3.56.5-2.5])结果如下ans=1.5000+1.6583i1.5000-1.6583i0.5000所以求得a=0.5，是正确的，而且在求解过程中只用到了y(n)与N下面来具体求解回声的这个问题利用自相关函数求N代码如下r=xcorr(y);plot(r);gridon;title('y的自相关函数');[u,v]=max(r);r1=r;r1(v-100:v+100,1)=0;[u1,v1]=max(r1);N=v-v1;下面是自相关函数的图像程序运行结果如下N=2400这与产生回声是所给的N是对的上的，求出的N是正确的，求得N后，就可以估计反射物的距离反射物的距离估计N=2400，而采样频率为fs=8000Hz，可得延时t=N/fs=0.3s设反射物距离为s，声速为v，则有2*s=v*t得s=v*t/2声速大概在340m/s左右，把t=0.3s,v=340m/s代入，求得s=51m，即反射物的距离为51米接着来求解a已知N，求解a上面已经说明了求解的方法，直接使用，代码如下fork=1:11t(k)=(-1)^k*y(2400*(k-1)+1,1);endroots(t)求解结果如下ans=4.8486+4.4069i4.8486-4.4069i-4.2880-0.5868+0.5978i-0.5868-0.5978i-0.60450.2255+0.6507i0.2255-0.6507i0.91810.5000把不合理的根舍掉，则只剩0.9181和0.5，这个时候就遇到了问题，有两个根都合理，该如何解决？在前面我们是从y(1)开始，一级级往后推，用的系数是y(1),y(2401),y(4801)...当然我们也可以从y(2)开始，一级级往后推，用系数y(2),y(2402),y(4802)...对于a来说是不变的，所以在两个求根结果中一样的根就是a，从y(2)开始的求解代码如下fork=1:11t(k)=(-1)^k*y(2400*(k-1)+2,1);endroots(t)求解结果如下ans=5.8416-5.2057-0.8448+0.3529i-0.8448-0.3529i0.1828+0.6540i0.1828-0.6540i0.63340.5000-0.1119与上面的结果比较，只有0.5是共同的根，所以a=0.5，这与在产生回声中给定的a是一样的，所以a的求解正确求得N与a后，下面就可以恢复原信号了从带有回声的声音信号中恢复原信号由上面的求解得到N=2400，a=0.5，所以y(n)=x(n)+0.4*x(n-2400)现在是已知y(n)，求解x(n)，这个时候x(n)是未知其系统函数为2400()1()()10.4*XzHzYzz知道系统函数后，可以调用filter函数filter是一维数字滤波器其使用方法如下：Y=filter(B,A,X)，输入X为滤波前序列，Y为滤波结果序列，B/A提供滤波器系数，B为分子，A为分母整个滤波过程是通过下面差分方程实现的：a(1)*y(n)=b(1)*x(n)+b(2)*x(n-1)+...+b(nb+1)*x(n-nb)-a(2)*y(n-1)-...-a(na+1)*y(n-na)下面从带有回声的声音信号中恢复原信号并保存，并且画出时域图与频域图代码如下a=[1,zeros(1,2399),0.5];b=[1];x1=filter(b,a,y);wavwrite(x1,'恢复后的信号')subplot(3,1,1);plot(x1);gridon;xlabel('时间');ylabel('幅值');title('恢复后的信号时域波形');subplot(3,1,2);wx1=fft(x1);f=(0:3*fs+N-1)*fs/(3*fs+N);plot(f,abs(wx1));gridon;xlabel('频率');ylabel('幅值');title('幅频特性');subplot(3,1,3);plot(f,angle(wx1));gridon;xlabel('频率');ylabel('相位');title('相频特性');恢复后的信号附件（双击打开）：恢复后的信号.wav图形如下：再与原始信号的图形（如下）进行比较几乎是一样的，这样就从带有回声的声音信号中恢复出了原信号说明：录制的声音音量好像较小，若查看语音文件，带上耳机可能会好些。现把三个语音文件附在下面（双击打开）：原始信号：原始信号.wav加回声后的信号：加回声后的信号.wav恢复后的信号：恢复后的信号.wav小结：这个过程大体上完成了所要求的功能：采集一个语音信号，加入回声，恢复原信号，估计反射物距离。在这里，有一个非常大的不足就是，对于各个函数都是直接引用已有函数，并未自己编程实现。在整个过程中，我认为有一点对于从回声信号中恢复原信号来说非常重要，那便是N与a的求解，而不应该把它们当做已知来直接滤波得到原序列。正因为要基于N与a的求解，带来了一些问题：如果a过小，则回声的加入对自相关函数来说影响不大，就很难求出N的值，所以带来了一定的局限性，这种基于N与a的求解的滤波相比较而言可能缺少一定的灵活性。参考文献：[1]程佩青编著.数字信号处理教程.北京：清华大学出版社,2007.[2]张德丰编著.详解MATLAB数字信号处理.北京：电子工业出版社，2010.[3]王彬主编，于丹汪洋副主编.MATLAB数字信号处理.北京：机械工业出版社，2010.