作者:师生共2人,其中,网安学院研究生朱可,教师张瑞
均衡算法有很多种,例如迫零均衡、MMSE均衡以及CMA均衡,它们各有优劣,可在不同情况下采用不同的算法。除此之外均衡还可分为时域均衡和频域均衡两种方式,时域均衡算法从冲激响应的角度进行改善,使其满足无码间串扰条件;而频域均衡算法从频率传递函数的角度进行改善,使其满足无失真传输条件。其中,迫零均衡和MMSE均衡算法既可以在时域实现,也可以在频域实现,而CMA均衡算法仅能在时域实现。
1 时域均衡算法
时域均衡算法的基本思想是将接收端接收到的信号通过滤波器进行处理,消除信道对信号的影响,从而还原原始信号。该算法通过对接收信号的每个符号进行处理,根据先前的符号推测当前符号,从而实现信道均衡。
时域均衡算法包括线性均衡算法和非线性均衡算法两种。
线性均衡算法通常采用卷积算法,通过反演信道响应的逆滤波器来抵消信道的影响,实现信道均衡。但是,由于卷积算法要求逆滤波器必须存在,因此该算法在实际应用中存在局限性。
非线性均衡算法则采用基于统计学习的方法,如最小均方误差(MMSE)等方法,通过学习信道的特征,对接收信号进行处理,消除信道对信号的影响。该算法可以适应更加复杂的信道环境,具有更好的性能和鲁棒性。
总的来说,时域均衡算法是数字通信系统中一种重要的信道均衡技术,可以有效抵消信道的影响,提高系统的可靠性和性能。时域均衡算法主要包括时域迫零均衡、时域MMSE均衡和CMA(Covariance Matrix Adaptation Evolution Strategy,协方差矩阵自适应演化策略)算法,以下分别对三种算法的基本原理进行说明。
1.1 时域迫零均衡
时域迫零均衡(Time-Domain Zero-Forcing Equalization Algorithm,TDE-ZF)的核心思想是通过调整均衡器的系数,使接收端的信号与发送端的信号在时域上尽量接近,从而实现信道均衡的目的。
如果接收端的位同步足够理想,那么采样将在最佳采样点位置,此时其他码元对当前采样点没有任何影响,当前时刻的取值仅由当前码元自身所决定,与其他码元无关,即无码间串扰。但是实际上位同步不会如此准确,而且每个码元经过信道后产生的波形也会发生一定的畸变,这使得当前采样点的值是前后多个码元波形在该点值的累加,即码间串扰。
时域迫零均衡消除码间串扰的原理是迫使前后几个码元对当前采样点的影响为零,具体操作就是将前后多个码元值进行加权累加,使其逼近当前时刻的理想值。其公式推导如下:
设发送的信号x(k)经过信道h(k)之后,受到噪声n(k)影响,则接收信号y(k)可以表示为:
线性均衡通常是利用滤波器实现的,对接收信号做均衡的操作是为了消除信道的影响,即将接收信号通过滤波器,得到能够判决的序列,其过程可以表示为:
其中c(k)表示均衡器的系数。
由于噪声无法预测,故忽略噪声,仅使用下式来趋近x(k),
上式的成立需要信道冲激响应和均衡器的系数在时域上满足以下条件:
那么对于时域迫零均衡只需要求出c(k)即可。
迫零均衡算法可以较为理想地消除码间串扰,而且均衡算法简单还有较低的复杂度,但是该算法没有考虑信道中存在的噪声。因此当某点衰落较严重时,均衡器会用较大增益对其补偿,但也会加大该点处的噪声。
1.2 时域MMSE均衡
在实际的通信过程中,信道往往是未知的,导致迫零均衡无法进行,而时域MMSE算法可以解决这个问题。
MMSE均衡算法的本质是让通过均衡器得到的信号值和真实值的均方差最小。设均衡后的信号值用S表示,信号的真实值用X表示,那么可以利用下式来建立基于MMSE均衡算法的延迟抽头滤波器:
其中表示均衡器的均衡系数矩阵,Y表示接收到的信号。
根据正交性原理,可以解出均衡系数矩阵,如下式所示:
其中,H表示信道的冲激响应,表示H的共轭转置矩阵,
和
分别表示发送信号X以及信道带外噪声N的平均功率,I表示的是单位矩阵。
当然,对于时域MMSE均衡而言,求出仍然无法直接计算出延迟抽头滤波器的加权值,此时需要作进一步推导,设其代价函数为下式:
其中,b(k)代表发送端信号也就是真实值,w代表均衡向量,r代表接收端前后几个信号组成的接收信号向量,对其进行求共轭导数并化简得到:
定义,
,R为接收信号的自相关矩阵,P为接收信号与发端信号的互相关矩阵。那么当导数为0时,可得:
由此得到时域MMSE均衡算法的均衡器。
MMSE算法不仅弥补了迫零均衡算法必须已知信道信息的劣势,拓展了均衡算法的应用范围,而且MMSE算法考虑到了噪声的影响,不会产生噪声增大的现象。
1.3 CMA算法
CMA算法实际上是经过三次改进后得到的结果,虽然MMSE均衡算法可以很好地实现理想中的均衡效果,但是MMSE均衡算法存在的三个缺点使其在实际中难以应用。
首先MMSE均衡算法需要对R矩阵进行求逆,而求逆操作在实际的FPGA中实现起来极为复杂,可通过最速下降法来解决这个问题;除此之外MMSE均衡算法只适应于慢衰落信道,MMSE均衡算法默认导频与数据经过的信道相同,因此无法适应快速变化的信道,在最速下降法基础上改进的LMS算法可以解决这个问题;最后,MMSE均衡算法需要使用导频头,需事先传输一段已知的序列来计算均衡器的抽头系数,这样会降低系统的传信率,而在LMS算法基础上改进的CMA算法可以解决这个问题。
对于最速下降法,其改进之处在于避免了求逆操作,通过迭代方式代替求逆,其核心思想是用不断迭代使均方误差下降,且迭代方向是均方误差下降最快的方向,由此得到最速下降法。
对二维曲面的某一点而言,其梯度一定是在该点上升最快的方向,那么负梯度值是在该点数值下降最快的方向,均方误差函数正是一个二元函数,其负梯度方向为其均方误差值下降最快的方向,首先得到梯度值,即:
其中,,
,b为发射端原始数据,r为接收端数据。然后假设已知当前时刻的均衡器w(n),希望下一时刻得到的均衡器w(n+1)是沿着下降最快的方向,走了一小步,其中
是步长。
将代价函数的梯度代入迭代公式,能得到:
按照上式进行计算,不需要求取R矩阵的逆变换,解决了MMSE算法的第一个难点。得到的迭代式子中只有步进量μ和初始滤波器w(0)是未知量,其他全部都是已知量,故接下来只需要确定μ和w(0)行 。
μ代表每一次迭代的步进长度,其值的选取需适当,否则过大会影响结果的收敛性,过小会导致收敛的时间过慢,故需根据实际情况选择合适的μ进行迭代。
对于w(0)初值的选取,由于对单边长度为L0的均衡器,其理想时域图中均衡器的最大值出现在第L0+1个点上,故在初值的选取上应取第L0+1位为1,其余位为零的序列。
对于LMS算法,其改进之处在于使算法适应快衰落信道,由于最速下降法是从最快方向迭代的,故其需要知道信道的统计特性,而LMS算法是在随机梯度方向迭代的,故其无需事先知晓信道统计特性,则即使信道快速变化对其也无太大影响。
将R和P的具体值代入最速下降法的公式得到下式:
对上式进行一定的改进,不再对p和R做统计平均,仅用当前时刻的值计算出随机梯度,来代替真实梯度,如下式所示:
对其进行化简得到下式:
其中y代表均衡后的数据,e代表原始数据与均衡后的数据的差值。
LMS算法和最速下降法相比,其能自适应信道变化,而且运算量小,但是由于其下降是随机的,故下降速度比最速下降法慢。
进一步改进去掉导频,通过传输信号幅值固定的特质来进行盲均衡的,即CMA算法。
导频本质是一段已知的期望函数,即b(n)。先用已知的导频将滤波器训练到收敛,然后滤波器可以被用于后面未知的信号。
首先考虑MMSE的标准形式解:
,其中,
,
,从上式可以看出,要想得到均衡器参数,必须首先知道导频b(n),而且对于MMSE均衡的代价函数而言,也需要真实信号来计算其均方误差,故对其代价函数进行一定修改,如下式所示:
其中R为模值,此算法适用于信号模值不变的通信系统,也可以经过一些转换应用到部分信号模值改变的通信系统中。对其进行求导并化简得到下式:
这里均衡器序列用c来表示,z为均衡后的值,结合梯度下降法的公式 可以得到CMA算法在梯度下降法中的迭代公式:
进一步,用随机梯度代替真实梯度,得到:
CMA算法收敛后,可以证明等于MMSE,其不需要b(n)导频,是盲估计。但是其收敛时间不能保证,主要依赖于信号的好坏。
2 频域均衡算法
频域均衡算法是一种数字信号处理技术,用于消除信号中的频率畸变或失真。在通信系统中,信号在传输过程中可能会受到信道衰落、多径效应等影响,导致信号在频域上出现失真。频域均衡算法通过对接收信号进行FFT(快速傅里叶变换)得到信号的频域表示,然后根据信道响应的估计值进行均衡,以消除信号的失真。
频域均衡算法与时域均衡算法相比,有一个显著优势,即频域均衡算法可以克服信道零点。时域均衡算法在多次仿真的过程中,总会出现均衡效果很差的情况,这时如果观察信道的幅频响应,会发现信道有一个接近零点的位置,如下图所示:
频域的集中等于时域的弥散,当频域有一个尖锐的深衰落,就意味着时域要接近无限长,但均衡器的阶数是有限的,因此这时,时域均衡往往不能起到很好的效果。而频域均衡的核心是在频域对信号进行均衡,不存在这种限制,故频域均衡算法比时域均衡算法更稳定,适用于多种信道条件。
2.1 频域迫零均衡
对于频域迫零均衡,其思想和时域迫零均衡是一致的,只是需要在频域考虑均衡,对于时域而言需要消除码间串扰,对应到频域是满足无失真传输条件。在上文时域迫零均衡中,推导出其在时域上需要的满足式,则在频域上需要满足:
那么均衡器频域表达式为:
这时仅需求出C(f)即可。
2.2 频域MMSE算法
频域MMSE算法的思想和时域是一致的,可以借此将其扩展到频域(已经过严格推导的证实),得到:
其中U为信道的频域特征,得到频域的均衡器后将其与FFT后的信号相乘,再取IFFT可以得到均衡后的信号。
3 代码参考
3.1 模拟信号发送和接收处理过程
% M为数据序列的block数量(单个信道状态下频域均衡考察的bit数量:M×128×2,时域约为M×160×2):
M = 400;
% 尝试多少种随机生成的信道(多径数、延迟和各径幅度恒定,各径复增益相角随机):
ChnStateTryNumber = 10;
% 时域均衡的阶数(单边):
L0 = 15;
% 导频序列:
bPilot = sqrt(2)*exp(1j*3*pi*(0:127).^2/128);
% 过采样率(复基带信号每个符号的采样点数):
SmpPerSymb = 8;
% 基带脉冲波形滚降系数:
Rolloff = 1;
htx = rcosdesign(Rolloff,6,SmpPerSymb,'sqrt');
% 循环前缀长度(按符号个数):
LengthOfCP = 32;
bPilotCP = bPilot(end-LengthOfCP+1:end);
% 加入循环前缀之后的导频序列:
bPilotCPextend = [bPilotCP,bPilot];
% 加入循环前缀之后的导频复基带信号:
sPilotBB = conv(upsample(bPilotCPextend,SmpPerSymb),htx);
CPG02ndpath = sqrt(4/5)*exp(1j*rand*2*pi);
CPG03rdpath = sqrt(2/3)*exp(1j*rand*2*pi);
CPG04thpath = sqrt(1/2)*exp(1j*rand*2*pi);
CPG05thpath = sqrt(2/5)*exp(1j*rand*2*pi);
CPG06thpath = sqrt(2/5)*exp(1j*rand*2*pi);
CPG07thpath = sqrt(1/3)*exp(1j*rand*2*pi);
CPG08thpath = sqrt(1/2)*exp(1j*rand*2*pi);
CPG09thpath = sqrt(1/4)*exp(1j*rand*2*pi);
CPG10thpath = sqrt(1/5)*exp(1j*rand*2*pi);
% % --- 预置多径模式(仅用于测试):
% load CPGdataNo02_SCTDE4Comparison_v1
% 经历多径衰落的复基带导频信号:
sPilotBBmultipath = ...
[sPilotBB,zeros(1,24)] + ... 主径(第一径)
[0,CPG02ndpath*sPilotBB,zeros(1,23)] + ... 第二径
[zeros(1,3),CPG03rdpath*sPilotBB,zeros(1,21)] + ... 第三径
[zeros(1,7),CPG04thpath*sPilotBB,zeros(1,17)] + ... 第四径
[zeros(1,9),CPG05thpath*sPilotBB,zeros(1,15)] + ... 第五径
[zeros(1,13),CPG06thpath*sPilotBB,zeros(1,11)] + ...第六径
[zeros(1,15),CPG07thpath*sPilotBB,zeros(1,9)] + ... 第七径
[zeros(1,18),CPG08thpath*sPilotBB,zeros(1,6)] + ... 第八径
[zeros(1,19),CPG09thpath*sPilotBB,zeros(1,5)] + ... 第九径
[zeros(1,24),CPG10thpath*sPilotBB];% 第十径
% 荷载数据的(payload)复基带QPSK信号经过与导频相同的衰落信道:
bTrafficImat = 1-2*randi([0 1],128,M);
bTrafficQmat = 1-2*randi([0 1],128,M);
bTrafficMat = bTrafficImat + 1j*bTrafficQmat;
bTrafficCPmat = bTrafficMat(end-LengthOfCP+1:end,:);
% 加入循环前缀之后的数据序列:
bTrafficCPextend = reshape([bTrafficCPmat;bTrafficMat],1,[]);
% 加入循环前缀之后的数据复基带信号:
sTrafficBB = conv(upsample(bTrafficCPextend,SmpPerSymb),htx);
% 多径衰落信道(假定数据经历的多径信道与导频完全一致):
% ※ 除0延迟主径外,另有9个从径,其延迟相对于主径见导频处定义。
sTrafficBBmultipath = ...
[sTrafficBB,zeros(1,24)] + ... 主径(第一径)
[0,CPG02ndpath*sTrafficBB,zeros(1,23)] + ... 第二径
[zeros(1,3),CPG03rdpath*sTrafficBB,zeros(1,21)] + ... 第三径
[zeros(1,7),CPG04thpath*sTrafficBB,zeros(1,17)] + ... 第四径
[zeros(1,9),CPG05thpath*sTrafficBB,zeros(1,15)] + ... 第五径
[zeros(1,13),CPG06thpath*sTrafficBB,zeros(1,11)] + ...第六径
[zeros(1,15),CPG07thpath*sTrafficBB,zeros(1,9)] + ... 第七径
[zeros(1,18),CPG08thpath*sTrafficBB,zeros(1,6)] + ... 第八径
[zeros(1,19),CPG09thpath*sTrafficBB,zeros(1,5)] + ... 第九径
[zeros(1,24),CPG10thpath*sTrafficBB];% 第十径
% =====================================================================
% 对导频的信号处理(无噪声,因此仿真可认为是在信道状态理想已知下进行的):
% =====================================================================
% 接收机匹配滤波:
rPilotBBpostMF = conv(sPilotBBmultipath,htx);
% 接收机去循环前缀,得到长度为1个block(128个符号)的观察窗口:
rPilotBBremoveCP = rPilotBBpostMF(151:1174);
% 位同步,找最佳采样点:
PeakFind = mean(reshape(abs(rPilotBBremoveCP),SmpPerSymb,[]).');
% 使用“眼图最大张开位置”作为最佳采样点:
PeakPosition = find(PeakFind==max(PeakFind));
% 抽取,得到位同步输出(1倍符号速率的最佳采样点序列):
SymbTimingDout = downsample(rPilotBBremoveCP(PeakPosition:end),SmpPerSymb);
% 频域信道估计:
ChnFDcoe = fft(SymbTimingDout)./fft(bPilot);
% =====================================================================
3.2 均衡算法代码
% 时域迫零均衡:
% =========================================================================
% 根据先验信道知识,生成信道矩阵,建立二进制符号序列和抽样判决序列的线性关系;
% 实际情况下,信道知识(“uMat”)需要靠信道估计(Channel Estimation)获得。
u = downsample(h(3:end),8);
u = u(end:-1:1);
uMat = [[u(3:end),0,0];[u(2:end),0];u;[0,u(1:end-1)];[0,0,u(1:end-2)];[0,0,0,u(1:end-3)]];
uMatInv = inv(uMat);
EqualizerVec = uMatInv(3,:);
% 结合卷积的定义思考一下:“EqualizerUIR”为什么等效于“EqualizerVec”的逆序:
EqualizerUIR = EqualizerVec(end:-1:1);
% =========================================================================
% 时域MMSE均衡:
% =====================================================================
% 位同步,找最佳采样点:
PeakFind = sum(reshape(abs(rTrafficBBpostMF(1:fix(length(rTrafficBBpostMF)/SmpPerSymb)*SmpPerSymb)),SmpPerSymb,[]).');
PeakFind = PeakFind - min(PeakFind);
% 使用“眼图最大张开位置”作为最佳采样点:
PeakPosition = find(PeakFind==max(PeakFind));
% 抽取,得到位同步输出(1倍符号速率的最佳采样点序列):
SymbTimingDout = downsample(rTrafficBBpostMF(PeakPosition:end),SmpPerSymb);
% 截掉头尾哪些幅度不够的数据(之所以前后各截掉7个,是多次观察后的经验所得):
SymbTimingDout = SymbTimingDout(8:end-8);
SymbTimingDout1 = SymbTimingDout(1:fix(length(SymbTimingDout)/L)*L);
rMatSize2 = length(SymbTimingDout)-2*L0;
rMatIdx = repmat((1:L).',1,rMatSize2) + repmat((0:rMatSize2-1),L,1);
rMat = SymbTimingDout(rMatIdx);
% 算数平均代替统计平均,求输入信号自相关矩阵R:
Rmat = zeros(L,L);
for q = 1:rMatSize2
Rmat = Rmat + rMat(:,q)*rMat(:,q)';
end
Rmat = Rmat/rMatSize2;
% 算数平均代替统计平均,求向量p:
% 期望响应:
d = bTrafficCPextend(L0+1:L0+rMatSize2);
dConjMat = repmat(conj(d),L,1);
p = mean((dConjMat.*rMat).').';
% 维纳解:w=inv(R)*p
wOpt = Rmat\p;
% 基于直接矩阵乘法的均衡:
MMSETDEdout = wOpt'*rMat;
%最速下降法
t=size(rMat,2);
ws=zeros(L,t+1);
ws(16,1)=1;
for n=1:t
ws(:,n+1)=ws(:,n)-0.01.*(Rmat*ws(:,n)-p);
end
ws_f=ws(:,n+1);
SDTDEdout = ws_f'*rMat;
%LMS
u=0.0002;
ws1=zeros(L,t+1);
ws1(16,1)=1;
for n=1:t-15
r=rMat(:,n);
b1=bTrafficCPextend(1,n+15);
b2=conj(b1);
ws1(:,n+1)=ws1(:,n)-u*r*(r'*ws1(:,n)-b2);
end
ws_f1=ws1(:,n+1);
%ws_f=ws_f(end:-1:1);
LMSTDEdout = ws_f1'*rMat;
zk=16;
dLMS = bTrafficCPextend(zk:end+zk-31);
% 基于CMA的自适应均衡:
% 转行向量:
EqualizerDin = SymbTimingDout;
% 对均衡器输入进行归一化幅度:
EqualizerDin = 1/sqrt(mean(abs(EqualizerDin).^2))*EqualizerDin;
muup = 0.001;
% 定义恒模算法的“模”:
R2 = 2;
cMat = zeros(L,length(EqualizerDin)-L+1);
% CMA均衡器向量赋初值:
cMat(:,1) = [zeros(L0,1);1;zeros(L0,1)];
% 迭代过程:
EqualizerDout = zeros(size(EqualizerDin)-L+1);
cMatSize2 = size(cMat,2);
CostFuncValue = zeros(1,cMatSize2);
for n = 1:cMatSize2
y_n = EqualizerDin(n:n+L-1).';
c_n = cMat(:,n);
z_n = y_n.'*c_n;
EqualizerDout(n) = z_n;
CostFuncValue(n) = (abs(z_n)^2-R2)^2;
c_nPlus1 = c_n - muup*conj(y_n)*z_n*(abs(z_n).^2-R2);
% 循环结束后,c_nPlus1即为稳态下的均衡器抽头系数:
if (n<cMatSize2)
cMat(:,n+1) = c_nPlus1;
end
end
EqualizerDoutSteady = conv(c_nPlus1(end:-1:1),EqualizerDin);
EqualizerDoutSteady = EqualizerDoutSteady(L:end-L+1);
% 频域迫零,使用导频辅助的信道估计结果进行频域均衡,再将均衡结果逆变换到时域:
EqualizerDoutZFmat = ifft(fft(SymbTimingDoutMat)./repmat(ChnFDcoe.',1,M));
EqualizerDoutZF = reshape(EqualizerDoutZFmat,1,[]);
% MMSE频域均衡:
% 求时域噪声功率:
NoisePastMF = conv(ChnNoise,htx);
NoisePastMFdownsampledMat = reshape(NoisePastMF(1:fix(length(NoisePastMF)/SmpPerSymb)*SmpPerSymb),SmpPerSymb,[]).';
NoisePastMFdownsampledPwr = var(NoisePastMFdownsampledMat);
SigmaSquared = mean(NoisePastMFdownsampledPwr);
% 生成MMSE频域均衡器:
ChnFDcoeConjugated = ChnFDcoe';
EqualizerMMSEvec = ChnFDcoeConjugated./(abs(ChnFDcoeConjugated).^2+SigmaSquared);
EqualizerDoutMMSEmat = ifft(fft(SymbTimingDoutMat).*repmat(EqualizerMMSEvec,1,M));
EqualizerDoutMMSE = reshape(EqualizerDoutMMSEmat,1,[]);
4 仿真结果
将七种均衡算法在QPSK调制,十径信道条件下进行仿真结果如下图所示。
从上图可见,在Eb/n0在1~14dB范围内,时域MMSE均衡、频域MMSE均衡、最速下降法、时域LMS均衡的误码率曲线几乎重合,且在Eb/n0为14dB时误码率可达到1e-5,体现了MMSE均衡方法的优越性,此外,最速下降法和时域LMS均衡可以较好代替时域MMSE均衡方法。
时域CMA算法、时域迫零均衡和频域迫零均衡算法性能接近,为达到相同误码率,所需Eb/n0比时域MMSE等算法高约0.8dB,即均衡性能差0.8dB。