均衡滤波算法
‘壹’ 自适应均衡器的原理是什么
均衡根据通信系统中的一项重要技术，分为两种方式：频域均衡和时域均衡。频域均衡是利用可调滤波器的频率特性来弥补实际信道的幅频特性和群延时特性，使包括均衡器在内的整个系统的总频率特性满足无码间干扰传输条件。时域均衡是直接从时间响应角度考虑，使包括均衡器在内的整个传输系统的冲激响应满足无码间干扰条件。频域均衡满足奈奎斯特整形定理的要求，仅在判决点满足无码间干扰的条件相对宽松一些。所以，在数字通信中一般采用时域均衡。
时域均衡器可以分两大类：线性均衡器和非线性均衡器。如果接收机中判决的结果经过反馈用于均衡器的参数调整，则为非线性均衡器；反之，则为线性均衡器。在线性均衡器中，最常用的均衡器结构是线性横向均衡器，它由若干个抽头延迟线组成，延时时间间隔等于码元间隔 。非线性均衡器的种类较多，包括判决反馈均衡器(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％以上，如果滤波器的阶数更大，则速度会提高更多。