(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，具体方法和原理如下：
加回声后的信号.wav
恢复后的信号：
恢复后的信号.wav
小结：
这个过程大体上完成了所要求的功能：采集一个语音信号，加入回声，恢复原信号，估计反 射物距离。在这里，有一个非常大的不足就是，对于各个函数都是直接引用已有函数，并未 自己编程实现。在整个过程中，我认为有一点对于从回声信号中恢复原信号来说非常重要，
H (z)
X (z) 1 Y ( z ) 1 0.4* z 2400
知道系统函数后，可以调用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); grid on; xlabel('时间'); ylabel('幅值'); title('恢复后的信号时域波形'); subplot(3,1,2); wx1=fft(x1); f=(0:3*fs+N-1)*fs/(3*fs+N); plot(f,abs(wx1)); grid on; xlabel('频率'); ylabel('幅值'); title('幅频特性'); subplot(3,1,3); plot(f,angle(wx1)); grid on; xlabel('频率'); ylabel('相位') ;title('相频特性');
数字信号处理课程设计
题目：基于 MATLAB 的回波信号的产生与消除 课程：MATLAatlab 采集一段语音，在这段语音的基础上，加入一定延时和衰 减的回音，最后消去回音并且测出延时时间来计算障碍物距离
正文
①设计目的与要求
采集语音：采集一段语音，绘制其时域波形，对此音频信号用 FFT 作谱分析。 加入回声：对采集的语音进行处理，加入一段回声，并绘制其时域波形，对其进行 FFT 频
利用自相关函数求 N
代码如下 r=xcorr(y);
plot(r); grid on; 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后，就可以估计反射物的距 离
加回声后信号附件(双击打开):
加回声后的信号.wav
(4) 加回声后信号的时域波形，FFT 频谱分析
代码如下： subplot(3,1,1); plot(y); grid on; xlabel('时间'); ylabel('幅值'); title('加回声后信号时域波形'); subplot(3,1,2); wy=fft(y); f=(0:3*fs+N-1)*fs/(3*fs+N); plot(f,abs(wy)); grid on; xlabel('频率'); ylabel('幅值'); title('幅频特性'); subplot(3,1,3); plot(f,angle(wy)); grid on; xlabel('频率'); ylabel('相位') ;title('相频特性'); 图如下：
现在逐级代入，可得 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)=[1 2 3 4 5] a=0.5 N=2 则有 x(n) =1 2 3 4 5 x(n-2)= 0.5 1 1.5 2 2.5 y(n) =1 2 3.5 5 6.5 2 2.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.5 6.5 -2.5]) 结果如下 ans = 1.5000 + 1.6583i 1.5000 - 1.6583i 0.5000 所以求得a=0.5，是正确的，而且在求解过程中只用到了y(n)与N 下面来具体求解回声的这个问题
如何求 N
利用自相关函数xcorr来估计N，对于信号x(n)，其长度为N，其求得的自相关函数为 r(m)=
x(n) * x(n m) ，其中m的范围为-(N-1)到N-1，而且显然是左右对称的。下面是一
n 1
N
个简单例子 运行代码 xcorr([1 2 3]) 结果如下 ans = 3.0000 8.0000 14.0000 8.0000 3.0000 自相关函数是对函数本身在两个时刻t1， t2的相关程度的一种衡量标准， 对于加回声后信号 y(n)=x(n)+ax(n-N),y(n)是由x(n)与它的一个衰减延时ax(n-N)叠加而成，因为相关函数是 函数本身相关程度的一种衡量，可以看到，y(n)的自相关函数将出现几个极值点，自相关函 数为r(m)=
0.1828 - 0.6540i 0.6334 0.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)是未知 其系统函数为
wx=fft(x); f=(0:3*fs-1)*fs/(3*fs); plot(f,abs(wx)); grid on; xlabel('频率'); ylabel('幅值'); title('幅频特性'); subplot(3,1,3); plot(f,angle(wx)); grid on; xlabel('频率'); ylabel('相位') ;title('相频特性'); 图如下：
恢复后的信号附件（双击打开）：
图形如下：
恢复后的信号.wav
再与原始信号的图形（如下）进行比较
几乎是一样的，这样就从带有回声的声音信号中恢复出了原信号
说明：
录制的声音音量好像较小，若查看语音文件，带上耳机可能会好些。现把三个语音文件附在 下面（双击打开）：
原始信号：
原始信号.wav
加回声后的信号：
把不合理的根舍掉，则只剩0.9181和0.5，这个时候就遇到了问题，有两个根都合理，该
如何解决？
在前面我们是从y(1)开始，一级级往后推，用的系数是y(1),y(2401),y(4801)... 当然我们也可以从y(2)开始， 一级级往后推， 用系数y(2),y(2402),y(4802)...对于a来说是 不变的，所以在两个求根结果中一样的根就是a，从y(2)开始的求解代码如下 for k=1:11 t(k)=(-1)^k*y(2400*(k-1)+2,1); end roots(t) 求解结果如下 ans = 5.8416 -5.2057 -0.8448 + 0.3529i -0.8448 - 0.3529i 0.1828 + 0.6540i
x(n) * x(n m) ，对于y(n)极值点应该出现在m=0，m=+-N，这时候相关程度