目標反射回波檢測算法及其FPGA實現 之一：算法概述
目標反射回波檢測算法及其FPGA實現之一：算法概述
前段時間，接觸了一個聲吶目標反射回波檢測的項目。聲吶接收機要實現的核心功能是在含有大量噪聲的反射回波中，識別出發射機發出的激勵信號的回波。我會分幾篇文章分享這個基於FPGA的回波識別算法的開發過程和原碼，歡迎大家不吝賜教。以下原創內容歡迎網友轉載，但請註明出處： https://www.cnblogs.com/helesheng。本文首先簡要介紹基於FPGA的互相關反射回波檢測算法的主要設計思路。
聲吶測距的原理非常簡單，如下圖所示，抹香鯨在水下發射特定的聲波信號，並監聽感興趣的對象（烏賊）反射的該聲波的回波。抹香鯨就可以根據回波的時間、幅度和相位的信息，“推算出”烏賊的方位、距離、速度和體表硬度等基本物理信息。
圖1 聲吶基本原理示意圖
這裏我們只討論最簡單的物理信息——距離。抹香鯨可以根據反射波到達的時間“計算”出和獵物之間的距離：L=V×Δt/2
圖2 聲吶測距示意圖
由上圖可知，聲吶測距的問題可以被簡化為：如何在接收到的信號中準確的定位出發射信號的反射回波。造成反射信號定位困難的原因有二：反射信號經過目標反射和長距離傳播後信號強度大大降低；大量的噪聲進一步降低了接收信號的信噪比。
下圖上方是對一段長度為64個采樣點的正弦信號（約6個周期）加blackman窗後得到的理想回波信號。下圖下方則是用matlab產生的仿真回波信號，其中加入了標準差為800白噪聲。可以看出，想不經過任何算法就得到反射回波的位置不是一件容易的事。
圖3 理論回波信號和加入了噪聲後的仿真回波信號
一種最直接的思路是對信號求短時傅裏葉變換（short-time FFT）。在Matlab中仿真這一過程的MATLAB代碼如下。
%%這個腳本用於驗證短時傅裏葉變換檢測回波的算法
i=1:60;
s=2000*sin(2*pi*i*6/60);%加了blackman窗的信號模板
s=[0,0,s,0,0];
s=s.*(blackman(64))‘;
sl=[zeros(1,100),s,zeros(1,100)];%兩邊擴展後的信號
figure;
subplot(2,1,1);
plot(sl);grid;
title(‘理論回波信號‘)
n=randn(1,264)*800;%方差為800的噪聲信號
slp=sl+n;%加入噪聲後的擴展信號
subplot(2,1,2);
plot(slp);grid;
title(‘加入白噪聲後的仿真回波實際信號‘)
[S,F,T,P]=spectrogram(slp,64,63,64,50E3);%slp為信號，64為加窗的長度，60為每次傅裏葉變換交疊部分的長度，64為傅裏葉變換的長度，50E3是采樣率。
figure;
surf(abs(S));
xlabel(‘時間‘)
ylabel(‘頻率‘)
zlabel(‘信號強度‘)
title(‘短時傅裏葉變換得到的時頻圖‘)
得到仿真回波信號的時頻圖如下所示。從圖中可以發現短時短時傅裏葉算法能夠成功的檢出回波出現的時間。 
圖4 短時傅裏葉變換得到的時頻圖
但采用短時傅裏葉變換搜索反射回波的辦法有以下缺點：
1、FFT算法的乘加計算的計算量為N*Log(N),例如，當窗口長度為64時，進行的乘加運算次數為64*6=384次MAC（乘加）運算。為了獲得最夠高的時域分辨率和目標定位精度，需要足夠高的短時傅裏葉窗口交疊率（overlap rate）。例如，為了獲得固定采樣率下的最高定位精度，每次采樣後都應該進行一次FFT，每個A/D采樣間隔要執行的計算量是384次MAC。
2、短時傅裏葉變換需要對信號加窗截斷為較短的長度分別進行FFT，以提高結果在時域的分辨率，而加窗會在截斷邊界處產生能量泄露效應。而這一點在上圖中就有所體現，在時頻圖的邊界處產生較大的幹擾，影響算法對回波位置的判斷。
3、短時傅裏葉變換只在激勵信號為正/余弦信號時才有效。如果激勵信號在頻域不是一根唯一的譜線（也就是非正/余弦信號），則反射回波中也會含有其他頻率分量，對其實施FFT就很難判斷哪些是激勵哪些是噪聲了。
一種更優的算法使通過計算激勵信號M(t)和反射回波s(t)的互相關：
(1)
來搜索反射回波在時間軸上的位置。其物理解釋是：在回波信號中搜索和激勵信號M(t)“最相似”的地方，從而定位哪裏最像回波信號。下圖是對沒有添加噪聲的理論回波信號相關的結果。
圖5 理論回波信號和激勵的互相關
顯然，由於激勵信號具有一定的周期性，導致互相關結果具有更加復雜的周期性結構。如果直接采用式（1）的結果來計算回波位置是不明智的，需要進一步通過積分環節來降低結果中的高頻周期撥動。當然，由於互相關結果R(τ)的直流分量為0，直接對其積分也不行。可以先對其求平方，以產生與R(τ)的功率成正比，且只有正值的信號，才能用於積分計算。我設計了如下所示的目標函數來表征回波信號與激勵信號的相似程度。
(2)
其中t0是最後一步積分求和窗口的長度，如果也取激勵信號的長度，即64，則得到如下圖所示的理想信號仿真結果。
圖6 對互相關結果求平方的仿真
圖7 對平方結果截斷積分的結果
采用和上述短時傅裏葉變換同樣嚴格的加性噪聲參數驗證這一算法的Matlab仿真代碼如下所示。
%這個腳本用於驗證回波卷積算法
i=1:60;
s=2000*sin(2*pi*i*6/60);%加了blackman窗的信號模板
s=[0,0,s,0,0];
s=s.*(blackman(64))‘;
plot(s);grid;
title(‘5K信號的64個點@50KSPS‘)
sl=[zeros(1,100),s,zeros(1,100)];%兩邊擴展後的信號
figure;
subplot(5,1,1);
plot(sl);grid;
title(‘理論回波信號‘)
n=randn(1,264)*800;%方差為800的噪聲信號
slp=sl+n;%加入噪聲後的擴展信號
subplot(5,1,2);
plot(slp);grid;
title(‘加入白噪聲後的仿真回波實際信號‘)
tp=conv(slp,s);%加入噪聲後的仿真實際信號和模板卷積（由於激勵信號左右對稱，這裏采用卷積函數conv()計算互相關）
subplot(5,1,3);
plot(tp);grid;
title(‘實際信號和信號模板卷積的結果‘)
tp1=tp.^2;%計算卷積結果的平方（能量）
subplot(5,1,4);
plot(tp1);grid;
title(‘卷積結果的平方‘)
w=ones(1,64);%矩形窗，用於通過卷積對歷史數據求和
tp2=conv(w,tp1);%歷史數據累加
subplot(5,1,5);
plot(tp2);grid;
title(‘歷史數據累加‘)
其運行結果如下圖所示。
圖8 對加入噪聲的回波信號采用互相關算法的仿真結果
可以看出最後一行的積分結果具有相當高的信噪比，能夠容易的分辨出加入噪聲（方差也是800）的回波信號中回波的位置。
另外，這種算法非常好的解決了短時傅裏葉變換中出現的問題：
1、相關算法的乘加計算量為N，大大小於FFT的計算量N*log(N)，當激勵信號長度為64個采樣點時，一個采樣周期中的計算量僅為64次MAC。當激勵信號為左右對稱的信號時，計算量可以進一步減少為N/2次乘加運算。
2、對激勵加窗函數後在計算相關可以有效的避免加窗帶來的泄露效應，從上圖的仿真結果中可以看出（2）式積分結果兩端不存在幹擾。
3、可對任何形式的激勵信號（如方波或鋸齒波）實時相關計算，其頻譜中即使存在多條譜線或連續譜也不會影響（2）式積分結果。
細心的讀者會發現所謂“短時傅裏葉方法”和這裏提出的“互相關搜索法”具有相通之處，用通俗的語言解釋：
傅裏葉變換的實質就是讓信號和不同頻率的正弦信號相關，以計算信號和它們的“相似度”，從而得到信號在不同頻率下的“功率密度”，並進一步得到頻譜的。而我們直接采用互相關法的優勢就在於，不在需要計算奈奎斯特頻率以內其他頻率正弦信號和信號的相關了。當然，由於我們搜索的就是回波中和激勵信號頻率相同的信號所在的位置，這種“互相關搜索法”自然能夠省去“短時傅裏葉方法”中大部分的計算量。
關於上述算法在FPGA中的實際實現，請關註後續博文“用Verilog-HDL摩爾型狀態機進行流程控制”。
關於A/D和D/A轉換即實驗平臺搭建的過程，請參看本系列的上一篇博文“互相關/卷積/FIR電路的實現”。
目標反射回波檢測算法及其FPGA實現 之一：算法概述