基于四二阶归一化统计量的盲均衡MATLAB仿真

clear all;
close all;
clc;

N=50000;
BitsTx = floor(rand(1,N)*2);%舍去小数部分
SymQAMtmp = reshape(BitsTx,4,N/4).';%重构矩阵,形成5000行4列
SymQAMtmptmp= bi2de(SymQAMtmp,'left-msb');%将二进制数转换成十进制数，等效于产生发送序列
QAMTable = [-3-3*j -3-j -3+3*j -3+j -1-3*j -1-j -1+3*j -1+j 3-3*j 3-j 3+3*j 3+j 1-3*j 1-j 1+3*j 1+j];%生成16QAM星座图
SymQAM = QAMTable(SymQAMtmptmp+1);%将发送序列映射到星座图上

figure(1);
plot(SymQAM,'.r');grid on;%输出发送序列的星座图
title('16QAM星座图');
axis([-6.0 6.0 -6.0 6.0]);

%Ch= [0.1 0.1 0.9 -0.1];%信道1
Ch= [0.05 -0.063 0.088 -0.126 -0.25 0.9074 0.25 0.126 0.038 0.088];%信道2

SymCh = filter(Ch,1,SymQAM);
%ChSym = awgn(SymCh,20,'measured');%加入高斯白噪声
ChSym =importdata('Chsym.mat');

figure(2);
plot(ChSym,'+r');%hold on;%打印均衡前信号，红色
title('均衡前信号');
axis([-6.0 6.0 -6.0 6.0]);grid on;

QQ=50;%设定估值的初值
X_WA_ABS_2 = mean((abs(ChSym(QQ:-1:1))).^2);
X_WA_2 = mean((ChSym(QQ:-1:1)).^2);
X_WA_ABS_4 = mean((abs(ChSym(QQ:-1:1))).^4);

rou = 0.999;%初始化设置
delta = 0.999;
lemda = 0.999;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u=0.001，第一次均衡 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0];%设置信道矩阵初值
u = 0.001;%设置步长
x_wa(14:-1:1) = ChSym(14:-1:1);%前14个数据直接取均衡器输入数据

for j=15:N/4
x_wa(j) = ChSym(j:-1:j-14)*conj(w');

s = conv(Ch,w);%求总的信道矩阵
isi(j) = (sum((abs(s)).^2)-(max(abs(s)))^2)/sum((abs(s)).^2);%求码间干扰
isi(j) = 10*log10(isi(j));%取dB单位

E_x_wa_4 = mean((x_wa(15:-1:1)).^4);
E_x_wa_2_2 = (mean((x_wa(15:-1:1)).^2))^2;
differ = E_x_wa_4/E_x_wa_2_2-3;

if(differ < 0)
sgn = -1;
else
sgn = 1;
end

w = w+conj(ChSym(j:-1:j-14))*u*sgn*mean((((abs(x_wa(j)))^2)*X_WA_ABS_2-X_WA_ABS_4+((abs(X_WA_2))^2))*x_wa(j)-X_WA_2*X_WA_ABS_2*conj(x_wa(j)))*((X_WA_ABS_2)^(-3));
X_WA_ABS_2 = lemda*X_WA_ABS_2+(1-lemda)*(abs(x_wa(j)))^2;
X_WA_2 = delta*X_WA_2+(1-delta)*(x_wa(j))^2;
X_WA_ABS_4 = rou*X_WA_ABS_4+(1-rou)*(abs(x_wa(j)))^4;
end%均衡算法部分，采用四、二阶归一算法

for k=1:4000%去除前4000个点，显示收敛后的点
x_wa(k) = 0;
end

figure(3);%输出u=0.001时，均衡后信号
plot(x_wa,'+r');grid on;
title('u=0.001时，均衡后信号');
axis([-6.0 6.0 -6.0 6.0]);

figure(4)
plot(isi);hold on;
title('不同步长u下ISI变化图');
xlabel('迭代次数');ylabel('ISI(dB)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u=0.005，第二次均衡 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0];%设置信道矩阵初值
u = 0.005;%设置步长
x_wa(14:-1:1) = ChSym(14:-1:1);

for j=15:N/4
x_wa(j) = ChSym(j:-1:j-14)*conj(w');

s = conv(Ch,w);%求总的信道矩阵
isi(j) = (sum((abs(s))

.^2)-(max(abs(s)))^2)/sum((abs(s)).^2);%求码间干扰
isi(j) = 10*log10(isi(j));%取dB单位

E_x_wa_4 = mean((x_wa(15:-1:1)).^4);
E_x_wa_2_2 = (mean((x_wa(15:-1:1)).^2))^2;
differ = E_x_wa_4/E_x_wa_2_2-3;

if(differ < 0)
sgn = -1;
else
sgn = 1;
end

%w = w+sgn*u*( ((abs(x_wa(j)))^2)*x_wa(j)*X_WA_ABS_2+x_wa(j)*((abs(X_WA_2))^2) )/((X_WA_ABS_2)^3)*conj(ChSym(j:-1:j-14))-sgn*u*( x_wa(j)*X_WA_ABS_4+conj(x_wa(j))*((X_WA_2)^2)*X_WA_ABS_2 )/((X_WA_ABS_2)^3)*conj(ChSym(j:-1:j-14));
w = w+conj(ChSym(j:-1:j-14))*u*sgn*mean((((abs(x_wa(j)))^2)*X_WA_ABS_2-X_WA_ABS_4+((abs(X_WA_2))^2))*x_wa(j)-X_WA_2*X_WA_ABS_2*conj(x_wa(j)))*((X_WA_ABS_2)^(-3));
X_WA_ABS_2 = lemda*X_WA_ABS_2+(1-lemda)*(abs(x_wa(j)))^2;
X_WA_2 = delta*X_WA_2+(1-delta)*(x_wa(j))^2;
X_WA_ABS_4 = rou*X_WA_ABS_4+(1-rou)*(abs(x_wa(j)))^4;
end%均衡算法部分，采用四、二阶归一算法

for k=1:4000%去除前4000个点，显示收敛后的点
x_wa(k) = 0;
end

figure(5)%输出u=0.005时，均衡后信号
plot(x_wa,'+r');grid on;
title('u=0.005时，均衡后信号');
axis([-6.0 6.0 -6.0 6.0]);

figure(4)
plot(isi,'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u=0.01，第三次均衡 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = [0 0 0 0 0 0 0 1 0 0 0 0 0 0 0];%设置信道矩阵初值
u = 0.01;%设置步长
x_wa(14:-1:1) = ChSym(14:-1:1);

for j=15:N/4
x_wa(j) = ChSym(j:-1:j-14)*conj(w');

s = conv(Ch,w);%求总的信道矩阵
isi(j) = (sum((abs(s)).^2)-(max(abs(s)))^2)/sum((abs(s)).^2);%求码间干扰
isi(j) = 10*log10(isi(j));%取dB单位

E_x_wa_4 = mean((x_wa(15:-1:1)).^4);
E_x_wa_2_2 = (mean((x_wa(15:-1:1)).^2))^2;
differ = E_x_wa_4/E_x_wa_2_2-3;

if(differ < 0)
sgn = -1;
else
sgn = 1;
end

%w = w+sgn*u*( ((abs(x_wa(j)))^2)*x_wa(j)*X_WA_ABS_2+x_wa(j)*((abs(X_WA_2))^2) )/((X_WA_ABS_2)^3)*conj(ChSym(j:-1:j-14))-sgn*u*( x_wa(j)*X_WA_ABS_4+conj(x_wa(j))*((X_WA_2)^2)*X_WA_ABS_2 )/((X_WA_ABS_2)^3)*conj(ChSym(j:-1:j-14));
w = w+conj(ChSym(j:-1:j-14))*u*sgn*mean((((abs(x_wa(j)))^2)*X_WA_ABS_2-X_WA_ABS_4+((abs(X_WA_2))^2))*x_wa(j)-X_WA_2*X_WA_ABS_2*conj(x_wa(j)))*((X_WA_ABS_2)^(-3));
X_WA_ABS_2 = lemda*X_WA_ABS_2+(1-lemda)*(abs(x_wa(j)))^2;
X_WA_2 = delta*X_WA_2+(1-delta)*(x_wa(j))^2;
X_WA_ABS_4 = rou*X_WA_ABS_4+(1-rou)*(abs(x_wa(j)))^4;
end%均衡算法部分，采用四、二阶归一算法

for k=1:4000%去除前4000个点，显示收敛后的点
x_wa(k) = 0;
end

figure(6)%输出u=0.01时，均衡后信号
plot(x_wa,'+r');grid on;
title('u=0.01时，均衡后信号');
axis([-6.0 6.0 -6.0 6.0]);

figure(4)
plot(isi,'g');
legend('u=0.001','u=0.005','u=0.01');