均衡濾波演算法
『壹』 自適應均衡器的原理是什麼
均衡根據通信系統中的一項重要技術，分為兩種方式：頻域均衡和時域均衡。頻域均衡是利用可調濾波器的頻率特性來彌補實際信道的幅頻特性和群延時特性，使包括均衡器在內的整個系統的總頻率特性滿足無碼間干擾傳輸條件。時域均衡是直接從時間響應角度考慮，使包括均衡器在內的整個傳輸系統的沖激響應滿足無碼間干擾條件。頻域均衡滿足奈奎斯特整形定理的要求，僅在判決點滿足無碼間干擾的條件相對寬松一些。所以，在數字通信中一般採用時域均衡。
時域均衡器可以分兩大類：線性均衡器和非線性均衡器。如果接收機中判決的結果經過反饋用於均衡器的參數調整，則為非線性均衡器；反之，則為線性均衡器。在線性均衡器中，最常用的均衡器結構是線性橫向均衡器，它由若干個抽頭延遲線組成，延時時間間隔等於碼元間隔 。非線性均衡器的種類較多，包括判決反饋均衡器(DFE)、最大似然(ML)符號檢測器和最大似然序列估計等。均衡器的結構可分為橫向和格型等。因為很多數字通信系統的信道(例如無線移動通信信道)特性是未知和時變的，要求接收端的均衡器必須具有自適應的能力。所以，均衡器可以採用自適應信號處理的相關演算法，以實現高性能的信道均衡，這類均衡器稱為自適應均衡器。
自適應均衡器的工作過程包含兩個階段，一是訓練過程，二是跟蹤過程。在訓練過程中，發送端向接收機發射一組已知的固定長度訓練序列，接收機根據訓練序列設定濾波器的參數，使檢測誤碼率最小。典型的訓練序列是偽隨機二進制信號或一個固定的波形信號序列，緊跟在訓練序列後面的是用戶消息碼元序列。接收機的自適應均衡器採用遞歸演算法估計信道特性，調整濾波器參數，補償信道特性失真，訓練序列的選擇應滿足接收機均衡器在最惡劣的信道條件下也能實現濾波器參數調整，所以，訓練序列結束後，均衡器參數基本接近最佳值，以保證用戶數據的接收，均衡器的訓練過程成功了，稱為均衡器的收斂。在接收用戶消息數據時，均衡器還不斷隨信道特性的變化連續地改變均衡器參數。
均衡器的收斂時間受均衡演算法、均衡器結構和信道特性的變化情況所決定。通常，均衡器需要通過重復性地周期訓練保證能夠一直有效地抑制碼間干擾。所以，用戶數據序列需要被分割成數據分組或時隙分段發送。
均衡器通常工作在接收機的基帶或中頻信號部分，基帶信號的復包絡含有信道帶寬信號的全部信息，所以，均衡器通常在基帶信號完成估計信道沖激響應和解調輸出信號中實現自適應演算法等
同時給你一個模擬例子：
自適應均衡器的模擬程序：
% Illustration of the conventional RLS algorithm
close all;
W=2.9;
Nexp=10;
N=2000;
Nmc=1; % Number of ensemble realizations
M=11;%抽頭系數
lambda=0.99;%遺忘因子
varv=0.001;%雜訊方差
h=zeros(3,1);%h的初始化
er=zeros(N,Nmc);%er的初始化
h(1)=0.5*(1+cos(2*pi*(1-2)/W));
h(2)=0.5*(1+cos(2*pi*(2-2)/W));
h(3)=0.5*(1+cos(2*pi*(3-2)/W));
% 學習曲線
hc=[0 h(1) h(2) h(3)]';
n0=7;
t=(1:N)';
for i=1:Nmc
y=sign(rand(N,1)-0.5);%輸入信號
v=sqrt(varv)*randn(N,1);%雜訊信號
x=filter(hc,1,y)+v;%信號混合
x=[zeros(M-1,1);x];%x矩陣
yd=zeros(N+M-1,1); %延遲信號初始化
e=yd;
yd(n0+M-1:N+M-1)=y(1:N-n0+1);
% CRLS 演算法
% Initialization
lambda=0.98;
P=(10^-3)*eye(M,M);
c=zeros(M,1);
g=c;
glambda=g;
% 迭代范圍
for n=M:M+N-1
xn=flipud(x(n-M+1:n));
glambda=P*xn;
alphal=lambda+conj(glambda')*xn;
g=glambda/lambda;a(n)=1-conj(g')*xn;
P=(P-g*conj(glambda'))/lambda;
P=(P+P')/2;
e(n)=yd(n)-conj(c')*xn;
c=c+g*conj(e(n));
end
end
eplot=e(M:M+N-1).^2;
subplot(2,1,1), plot(t,abs(eplot))
ylabel('|e(n)|^2');
xlabel('n');
subplot(2,1,2), plot(t,a(M:M+N-1));
ylabel('\alpha(n)');
xlabel('n');
『貳』 減少雜訊的匹配濾波演算法
（1）傳統匹配濾波演算法
Rickett et al.（2001）給出了匹配濾波簡要的公式及運算元長度設計標准，本節給出了更為詳細的匹配 濾波公式，並給出推導公式基本條件和結果。
設同一地區不同時期Y1，Y2得到的地震數據分別為GY1（t），GY2（t），取Y1年份的地震記錄為參
考地震道，使Y2年份相應的地震記錄與之匹配。選取歸一化運算元p使得目標泛函：
海上時移地震油藏監測技術
極小。最終得到關於求解匹配濾波器｛P（m），m＝1，2，…，L｝的L個方程的方程組：
海上時移地震油藏監測技術
為意義更明確，對上面的公式進一步簡化，令
海上時移地震油藏監測技術
上兩式中：RY2Y2（m－n）為時間延遲為m－n的時期Y2地震記錄在設計窗口中的自相關；RY1Y2（n）為時間延遲為n的時期Y1與時期Y2地震記錄在設計窗口中的互相關，於是方程（4.8）可以進一步寫成：
海上時移地震油藏監測技術
求解方程組（4.11）得到匹配濾波器運算元｛P（m），m＝1，2，…，L｝，用
海上時移地震油藏監測技術
校正相應的地震剖面。通過實際數據處理結果驗證了上述推導的正確性和方法的有效性。
方程（4.11）寫成矩陣形式：
海上時移地震油藏監測技術
式中：M為時期Y2地震記錄在設計窗口中的自相關序列組成的Toeplitz矩陣，R為時期Y1與時期Y2地 震記錄在設計窗口中的互相關序列向量。求解方程（4.13）可採用Levinson遞推演算法，計算效率高。
為了減少噪音的影響，通常引入阻尼項，方程（4.13）變為
海上時移地震油藏監測技術
式中：μ為很小的數，通常為可設為0.01或0.001。
實際應用中，可以發現式（4.13）受雜訊的影響很大，不穩定。雖然加入阻尼項後結果有所改善，但 如何選取合適阻尼因子又是一個難題。為此推導新的匹配濾波表達形式，尋求更穩健的求解方法。
（2）新匹配濾波公式
同樣設同一地區不同時期Y1，Y2得到的地震數據分別為GY1（t），GY2（t），取Y1年份的地震記錄 為參考地震道，使Y2年份相應的地震記錄與之匹配。則匹配過程可描述為
海上時移地震油藏監測技術
其中M為GY2組成的褶積矩陣。如果設地震道的采樣點數為n，設計濾波器f長度為m，M則為（2×n－1）×m矩陣，為保持矩陣維數相同，一種方法是將GY1後面補零為（2×n－1）×1向量，另一種方法是取 矩陣M的前n×m項。如果採用第一種方法，可以驗證得到的公式與（4.13）式相同。在此採用後一種方 法，得到新的匹配濾波方程。只要設計濾波器f足夠長，總能滿足能量差e（f）最小，根據范數定義：
海上時移地震油藏監測技術
求解能量差e（f）最小問題可轉化為
海上時移地震油藏監測技術
即對濾波因子向量求導，最終可歸結為求解線性方程：
海上時移地震油藏監測技術
如果記A＝MTM，b＝MTGY1，方程（4.18）轉化為
海上時移地震油藏監測技術
（4.19）式形式上與（4.13）式類似，內容不同，不再是Toeplitz矩陣，因此不能應用Levinson遞推演算法求解。因此，引入奇異值分解方法求解方程（4.19）。
（3）基於奇異值分解的匹配濾波演算法
矩陣的奇異值分解，是矩陣計算中一套很有用的技術。它可以有效地處理系數矩陣是奇異的或者接 近奇異的方程組。對於矩陣A，如果A∈Rm×n，並且A的秩為r，總有
海上時移地震油藏監測技術
其中， V為正交陣。 ，並且 為A 的奇異值。
公式（4.20）即為矩陣A的奇異值分解，根據正交矩陣的性質：
海上時移地震油藏監測技術
很容易表示出矩陣A的逆矩陣
海上時移地震油藏監測技術
將式（4.22）帶入式（4.19）中，得到濾波因子的表達式為
海上時移地震油藏監測技術
實際計算中，當A是奇異陣出現奇異值，或A接近奇異或病態矩陣時，（4.23）式的計算過程就無法進行。這時可將出現的奇異項 （σk是零，或者數值很小）簡單地替換成零或很小的常數，通過這種方法能得 到方程穩定的解。
對於實際含有雜訊的信號，信號能量主要分布在奇異值大的分量上，因此去除小奇異值同時能消除 雜訊影響。通常可選取某一能量百分比的奇異值作為去除的閾值，以這種方式既能克服A接近奇異或病 態矩陣的影響，又能減小雜訊的影響，使濾波因子穩健。
（4）模擬數據驗證
模擬得到一組存在時間、振幅、頻率、相位差異的信號，作為基測線與監測測線地震道，對監測測 線地震道加入不同比例的隨機雜訊，組成驗正演算法有效性的數據體，如圖4.10所示。分別用傳統的匹配 濾波方法和重新推導的基於奇異值分解的匹配濾波方法進行匹配處理，比較匹配後基測線與監測測線振 幅差異，結果見圖4.11和圖4.12。可以看出，傳統匹配濾波公式的計算結果受雜訊的影響很大，而基於 奇異值分解的匹配濾波方法具有很好的抗雜訊能力。
圖4.10 模擬地震記錄（從上至下依次為加入0％，10％，20％，30％雜訊的信號）
圖4.11 傳統方法匹配結果
圖4.12 基於奇異值分解方法匹配結果
（5）實際數據驗證
選擇一塊同一地區兩次不同時間測得的兩條二維測線；選取油藏上方時間長度為300ms的窗口作為 濾波因子設計窗口，並以抽取其中139道構成驗證互均衡演算法的數據體（圖4.13，圖4.14）。分別採用 傳統匹配濾波公式與基於奇異值分解的匹配濾波兩種方法進行校正。比較差異剖面的平均能量，結果見 圖4.15。從圖中可知基於奇異值分解的匹配濾波方法具有更好的抗雜訊能力，匹配誤差遠小於傳統匹配 濾波。
圖4.13 某地區時間1地震記錄
圖4.14 某地區時間2地震記錄
圖4.15 兩種匹配方法結果誤差能量對比圖
本節推導了新的匹配濾波方程，提出基於奇異值分解的匹配濾波演算法，理論和實際數據都驗證了該 方法有效性。這里從計算精度上比較兩種匹配濾波演算法，實際處理時移地震數據時還要考慮計算時間，此時尋求快速的奇異值分解演算法是一種提高處理效率的方式，另外針對不同信噪比，將傳統匹配濾波算 法與基於奇異值分解的匹配濾波演算法結合應用同樣是一種很好的方式。總之，基於奇異值分解的匹配濾 波提高了匹配精度，有利於為時移地震解釋提供一致性更好的地震資料。
『叄』 LMS自適應演算法分析及在數字濾波器設計中的應用
自適應過程一般採用典型LMS自適應演算法，但當濾波器的輸入信號為有色隨機過程時，特別是當輸入信號為高度相關時，這種演算法收斂速度要下降許多，這主要是因為輸入信號的自相關矩陣特徵值的分散程度加劇將導致演算法收斂性能的惡化和穩態誤差的增大。此時若採用變換域演算法可以增加演算法收斂速度。變換域演算法的基本思想是：先對輸入信號進行一次正交變換以去除或衰減其相關性，然後將變換後的信號加到自適應濾波器以實現濾波處理，從而改善相關矩陣的條件數。因為離散傅立葉變換�DFT本身具有近似正交性，加之有FFT快速演算法，故頻域分塊LMS�FBLMS演算法被廣泛應用。
FBLMS演算法本質上是以頻域來實現時域分塊LMS演算法的，即將時域數據分組構成N個點的數據塊，且在每塊上濾波權系數保持不變。其原理框圖如圖2所示。FBLMS演算法在頻域內可以用數字信號處理中的重疊保留法來實現，其計算量比時域法大為減少，也可以用重疊相加法來計算，但這種演算法比重疊保留法需要較大的計算量。塊數據的任何重疊比例都是可行的，但以50％的重疊計算效率為最高。對FBLMS演算法和典型LMS演算法的運算量做了比較，並從理論上討論了兩個演算法中乘法部分的運算量。本文從實際工程出發，詳細分析了兩個演算法中乘法和加法的總運算量，其結果為：
復雜度之比＝FBLMS實數乘加次數／LMS實數乘加次數＝（25Nlog2N＋2N－4）／[2N(2N－1)]�
採用ADSP的C語言來實現FBLMS演算法的程序如下：
for(i＝0;i＜＝30;i＋＋)
{for(j＝0;j＜＝n－1;j＋＋)
{in[j]＝input[i×N＋j;]
rfft(in,tin,nf,wfft,wst,n);
rfft(w,tw,wf,wfft,wst,n);
cvecvmlt(inf,wf,inw,n);
ifft(inw,t,O,wfft,wst,n);
for(j＝0,j＜＝N－1;j＋＋)
{y[i×N＋j]＝O[N＋j]．re;
e[i×N＋j]＝refere[i×N＋j]－y[i×N＋j];
temp[N＋j]＝e[i×N＋j;}
rfft(temp,t,E,wfft,wst,n);
for(j＝0;j＜＝n－1;j＋＋)
{inf_conj[j]＝conjf(inf[j]);}��
cvecvmlt(E,inf_conj,Ein,n);
ifft(Ein,t,Ein,wfft,wst,n);
for(j＝0;j＜＝N－1;j＋＋)
{OO[j]＝Ein[j]．re;
w[j]＝w[j]＋2＊u＊OO[j];}��
}
在EZ－KIT測試板中，筆者用匯編語言和C語言程序分別測試了典型LMS演算法的運行速度，並與FBLMS演算法的C語言運行速度進行了比較，表2所列是其比較結果，從表2可以看出濾波器階數為64時，即使是用C語言編寫的FBLMS演算法也比用匯編編寫的LMS演算法速度快20％以上，如果濾波器的階數更大，則速度會提高更多。