1、信號(hào)與系統(tǒng)A（2）課程項(xiàng)目實(shí)施報(bào)告 題目( B )： 回聲信號(hào)檢測(cè)系統(tǒng)設(shè)計(jì) 組 號(hào)： B28 任課教師： 組 長：_ _成 員：_成 員： _成 員： _成 員：_聯(lián)系方式： 二零 15年 5 月 22 日目錄：一、課程項(xiàng)目實(shí)施方案二、項(xiàng)目的設(shè)計(jì)過程單次回聲的分析與消除1、回聲的產(chǎn)生及信號(hào)的時(shí)域、頻域分析2、自相關(guān)函數(shù)分析回聲延時(shí)量3、用迭代法求回聲衰減系數(shù)4、回聲的消除多次回聲的分析，與消除1、回聲的產(chǎn)生和回聲延時(shí)量的分析2、用歸一化法估算回聲衰減系數(shù)3、回聲的消除三、項(xiàng)目的結(jié)果分析四、分工合作情況參考文獻(xiàn)附錄（matlab代碼）一、 課程項(xiàng)目實(shí)施方案題目要求：B組

2、 回聲信號(hào)檢測(cè)系統(tǒng)設(shè)計(jì)設(shè)計(jì)一個(gè)回聲信號(hào)檢測(cè)系統(tǒng)，能夠檢測(cè)出原信號(hào)中包含的回聲信號(hào)，并分析回聲信號(hào)的參數(shù)。要求：（1）輸入混有回聲的聲音（或信號(hào)源）；（2）設(shè)計(jì)檢測(cè)系統(tǒng)，分析回聲信號(hào)產(chǎn)生的原因、分析回聲信號(hào)的幅度、分析回聲信號(hào)的延遲量； (3) 能夠現(xiàn)場(chǎng)演示。設(shè)計(jì)思想：我們組的設(shè)計(jì)思想是利用matlab軟件平臺(tái)，用程序人為制造回聲，并且在僅知道回聲信號(hào)的條件下，分析延時(shí)量和衰減系數(shù)，從而消除回聲，并且測(cè)出聲源與障礙物的距離。二、 項(xiàng)目的設(shè)計(jì)過程單次回聲的分析與消除1、回聲的產(chǎn)生及信號(hào)的時(shí)域、頻域分析對(duì)于已經(jīng)采集好的一段語音信號(hào)，利用matlab的”audioread”命令可以對(duì)該信號(hào)進(jìn)行讀取和采

3、樣，本項(xiàng)目中采用的文件為wav文件，采樣頻率為44100.x,fs=audioread(E:ceshi.wav);利用plot和fft函數(shù)可以對(duì)已采集的信號(hào)分別繪制時(shí)域波形圖和經(jīng)快速傅里葉變換得到的頻譜圖；abs和angle函數(shù)可以繪制信號(hào)頻域的幅值與相位圖。 利用zeros函數(shù)創(chuàng)建空矩陣可以對(duì)原始信號(hào)進(jìn)行延時(shí)，并產(chǎn)生回聲。本例子中，假設(shè)延時(shí)量為8000個(gè)采樣點(diǎn)，衰減系數(shù)為0.4。同樣對(duì)回聲信號(hào)y進(jìn)行時(shí)域、頻域分析。2、自相關(guān)函數(shù)分析回聲延時(shí)量自相關(guān)函數(shù)的數(shù)學(xué)定義為：同一時(shí)間函數(shù)在瞬時(shí)t和t+a的兩個(gè)值相乘積的平均值作為延遲時(shí)間t的函數(shù)，它是信號(hào)與延遲后信號(hào)之間相似性的度量。簡(jiǎn)單地說，自相關(guān)函

4、數(shù)是用來表征一個(gè)信號(hào)在不同時(shí)刻t1,t2狀態(tài)之間的相關(guān)程度，其在某一時(shí)間點(diǎn)的值越大，相關(guān)程度越高。在matlab中，對(duì)于一個(gè)混有回聲的語音信號(hào)，可以利用xcorr自相關(guān)函數(shù)可以求出回聲信號(hào)與原始信號(hào)的采樣點(diǎn)之差，進(jìn)而求得延時(shí)量。 從自相關(guān)函數(shù)圖像中可以看到，在采樣點(diǎn)的中心有一最大的幅值，其含義就是當(dāng)延時(shí)量為零時(shí)，信號(hào)相關(guān)程度最高；而在這一峰值的兩側(cè)有兩個(gè)次大的幅值，其含義就是在該處信號(hào)相關(guān)程度次高；而在別的點(diǎn)處其相關(guān)程度相對(duì)于該2個(gè)峰值都可以忽略，因此算出這2個(gè)峰值所對(duì)應(yīng)的采樣點(diǎn)之差，就是回聲信號(hào)中，原始信號(hào)與回聲信號(hào)對(duì)應(yīng)的采樣點(diǎn)之差。由此，可以計(jì)算出該信號(hào)中采樣點(diǎn)之差為8000個(gè)單位，這與我

5、們所假設(shè)的值是吻合的。最后，將這個(gè)采樣點(diǎn)之差除以采樣頻率，即，便可得到實(shí)際的延時(shí)量。3、用迭代法求回聲衰減系數(shù)已知回聲的延遲量后，還需要知道回聲的衰減系數(shù)才可以確定回聲信號(hào)。對(duì)于一個(gè)一次回聲信號(hào)，我們可以采用迭代的方法來求衰減系數(shù):設(shè)xn為原始信號(hào)，yn為原始信號(hào)與回聲信號(hào)的疊加則有: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序列的第一個(gè)點(diǎn)是沒有回聲的，在第（1+N）、（1+2N

6、）、（1+（k-1）N）的點(diǎn)都是由原始信號(hào)與衰減后回聲信號(hào)疊加而成的，最后一個(gè)（1+kN）個(gè)點(diǎn)是僅含衰減后的回聲的。由圖像可以看出:設(shè)y(n)的長度為L，對(duì)于k,則滿足1+kNL（且k是整數(shù)）這個(gè)k是唯一的并且可求的。因此由最后一個(gè)式子逐級(jí)帶回上面的式子，將含x(n)的項(xiàng)全迭代成含y(n)和衰減系數(shù)的項(xiàng)，最后可以得到一個(gè)關(guā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函數(shù)可以對(duì)進(jìn)行求解。但是對(duì)于求解出的k個(gè)根，排除虛數(shù)和不在01范圍內(nèi)的數(shù)，可

7、能仍有一些根是滿足條件的。這時(shí)候我們可以將上述的迭代公式從y(2)開始計(jì)算，即：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)=ax(2+(k-1)N)對(duì)于同一個(gè)信號(hào)衰減系數(shù)必定相同，那么在兩個(gè)方程中，相同的根便是所求的。我們組為了實(shí)現(xiàn)程序的的普遍性，對(duì)程序進(jìn)行了優(yōu)化，自動(dòng)找出兩組解中相同的解。因此用迭代法可以精準(zhǔn)的算出回聲的衰減系數(shù)。4、回聲的消除由上面求解得到的延時(shí)量N_delay=8000，衰減系數(shù)=0.4，可以寫出原始信號(hào)與含回聲信號(hào)的關(guān)系式y(tǒng)(

8、n)=x(n)+0.4x(n-8000),其系統(tǒng)函數(shù)為H(z)=在matlab中可以用filter函數(shù)對(duì)y信號(hào)進(jìn)行濾波：z即為濾波完后的信號(hào)，sound(z,fs)發(fā)現(xiàn)單次回聲已經(jīng)全部消除。多次回聲的分析，與消除1、回聲的產(chǎn)生和回聲延時(shí)量的分析多次回聲的產(chǎn)生原理與一次回聲類似，利用zeros函數(shù)創(chuàng)建空矩陣可以對(duì)原始信號(hào)進(jìn)行延時(shí)。本例中延時(shí)量分別為8000個(gè)和12000個(gè)采樣點(diǎn)，衰減系數(shù)分別為0.4和0.2.同樣的，用xcorr自相關(guān)函數(shù)計(jì)算延遲量。 其中N_delay1是第一次回聲延遲的采樣點(diǎn)，N_delay2是第二次回聲延遲的采樣點(diǎn)。2、用歸一化法估算回聲衰減系數(shù)加入兩次回聲后，無法采用之前

9、單次回聲的方法計(jì)算衰減系數(shù)，因?yàn)闊o法確定唯一的k來列出迭代的方程式。我們采用自相關(guān)函數(shù)的歸一化法估算回聲衰減系數(shù)。具體過程如下：設(shè)Ry(n)為y(n)序列的自相關(guān)函數(shù)，Rx(n)為x(n)序列的自相關(guān)函數(shù)。由自相關(guān)函數(shù)的定義可以化簡(jiǎn)得到以下方程：Ry(n)=(1+)Rx(n) +Rx(n-) +Rx(n+)Rx(n-) +Rx(n+)+Rx(n+-) +Rx(n+)將n=0,N2分別代入R-1(0)=(1+)Rx(0)+Rx(-) +Rx()Rx(-) +Rx()+Rx(-) +Rx() （2）考慮到Rx(-) =Rx()約等于0，Rx(-) =Rx()約等于0，Rx(-) =Rx()約等于0

10、，（2）式化為Ry (0)=(1+)Rx(0) （3）同理有Ry()約等于Rx(0) ，Ry()約等于Rx(0)因此有, 還可推廣到n次回聲由此可以估算出第K次衰減系數(shù)。 利用matlab的fsolve解非線性方程組的功能，便可估算出2個(gè)回聲的衰減量。與理論的0.4和0.2有5%左右的誤差。3、回聲的消除同樣利用filter函數(shù)。H(z)=執(zhí)行sound(z,fs)后，回聲明顯不見了，可以達(dá)到效果。三、 項(xiàng)目的結(jié)果分析利用回聲測(cè)距上述兩例單次回聲和多次回聲的消除效果是良好的，消除后的聲音完全聽不出回聲。此外通過對(duì)回聲的分析可以應(yīng)用與測(cè)距。其原理很簡(jiǎn)單：通過自相關(guān)函數(shù)我們已求出回聲的延遲采樣點(diǎn)個(gè)

11、數(shù)，N_delay。其真實(shí)時(shí)間其中Fs是采樣率，t的單位為秒。則聲源與障礙物的距離。以報(bào)告中的例子舉例：?jiǎn)未位芈曋蠳_delay=8000，F(xiàn)s=44100,則可算出t=0.18s，那么聲源與障礙物的距離s=340*0.18/2=30.6m四、分工合作情況查找資料：盛亦池 胡斌 劉楊雨編寫代碼：盛亦池 應(yīng)祺超代碼調(diào)試：胡斌 應(yīng)祺超 馬嘉敏終期報(bào)告：馬嘉敏 劉楊雨參考文獻(xiàn)：1、自相關(guān)函數(shù) 百度百科2、matlab課程設(shè)計(jì)-基于MATLAB的回波信號(hào)的產(chǎn)生與消除 道客巴巴3、信號(hào)與系統(tǒng)實(shí)驗(yàn)指導(dǎo)4、基于MATLAB的數(shù)字信號(hào)處理與實(shí)踐開發(fā)，王佳梅，西安電子科技大學(xué)出版社附錄一次回聲matlab代碼：

12、%制造回聲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;%原始信號(hào)時(shí)域、頻域、幅值、相位圖n1=0:22050;figure(1)subplot(2,2,1);plot(x);title(原始信號(hào)時(shí)域圖);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/

13、4),X(1:44100/4);title(原始信號(hào)頻譜圖);subplot(2,2,3);plot(abs(X);title(原始信號(hào)幅值);subplot(2,2,4);plot(angle(X);title(原始信號(hào)相位);%疊加信號(hào)的時(shí)域、頻域、幅值、相位圖figure(2);subplot(2,2,1);plot(y);title(帶有回聲回聲時(shí)域圖);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:

14、44100/4);title(帶有回聲信號(hào)頻譜圖);subplot(2,2,3);plot(abs(Y);title(回聲幅值);subplot(2,2,4);plot(angle(Y);title(回聲相位);%自相關(guān)函數(shù) 算延時(shí)量cy=xcorr(y);figure(3)plot(cy);title(多次回聲的自相關(guān)函數(shù)圖像);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:

15、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; endendfor 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)=

16、;jie2(find(jie2(:,1)=0)= ;jie1=jie1;jie2=jie2;jie=(intersect(int16(jie1*1000),int16(jie2*1000);format shortfinaljie=double(jie)/1000;%消除回聲a=zeros(1,N);a(1)=1;a(N_delay+1)=finaljie;z=filter(1,a,y);多次回聲matlab代碼：x,fs=audioread(E:ceshi.wav);%制造回聲N=length(x);x1=x(1:N);x2=x(1:N);x3=x(1:N);x1=x1,zeros(1,12

17、000);x2=zeros(1,8000),0.4*x2,zeros(1,4000);x3=zeros(1,12000),0.2*x3;y=x1+x2+x3;%原始信號(hào)時(shí)域、頻域、幅值、相位圖n1=0:22050;figure(1)subplot(2,2,1);plot(x);title(原始信號(hào)時(shí)域圖);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(原始信號(hào)頻譜圖);subplot

18、(2,2,3);plot(abs(X);title(原始信號(hào)幅值);subplot(2,2,4);plot(angle(X);title(原始相位);%疊加信號(hào)的時(shí)域、頻域、幅值、相位圖figure(2);subplot(2,2,1);plot(y);title(帶有回聲信號(hào)時(shí)域圖);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(帶有回聲信號(hào)頻譜圖);subplot(2,2,3);plot(abs(Y);title(幅值);subplot(2,2,4);plot(angle(Y);title(相位);%自相關(guān)函數(shù)，計(jì)算延時(shí)量cy=xcorr(y);figure(3)plot(cy);title(多次回聲的自相關(guān)函數(shù)圖像);h1_y,h1=max(cy);cy2=cy;cy2(1,h1-1000:h1+1000)