cma.m

描述恒模算法的一个仿真程序

function []=cma(Le)
%本例计算均衡器长度与误码率的关系，信道是chan7的一段，h=C(20:40);
sta=randn('state');
clc;

P=1;%oversample ratio
sgma=0.0562%.01;%.05;%0.001;%0.0562;%相当于20dB
Lsh=6;%设定的信道长度
N=13;%N,表示累计的接收数据量;
%N只是在接收时起作用
Nu=5000;%发送的符号数目
Neu=4000;%CMA算法用到的数据长度
N0=50*P;%开始的符号序列数，从20个符号开始恢复；
if P==1;
    N0=30;
end
% hl=getchanel2(0.45,1.9,Lsh)+0.42*getchanel2(0.45,3.05,Lsh)+0.08*getchanel2(0.45,5.3,Lsh);

hl=getchanel2(0.45,1.9,Lsh)+0.42*getchanel2(0.45,3.05,Lsh)+0.53*getchanel2(0.5,5.3,Lsh);
hl=hl/norm(hl,'inf');
hz=hl(1:100:end);%没有过采样时的离散信道；
Lh=length(hz)-2;
ht=hl(1:100/P:end);%过采样率为L的离散信道；
h=reshape(ht(1:end-1),P,Lsh);

%-----------------------------------------------------------------------------------
Zero=zeros(P,1);Hc=[h zeros(P,N-1)];TR=diag(ones(1,Lh+N-1),1);TR(Lh+N,1)=1;H=Hc;
for i=2:N;%形成式2.7中的toeplitz矩阵H
    Hctemp=Hc*TR^(i-1);
    H=cat(1,H,Hctemp);
end
% load unablesolve;


hl=hl/norm(hl,'inf');
hz=hl(1:100:end);%没有过采样时的离散信道；
Lh=length(hz)-2;
ht=hl(1:100/P:end);%过采样率为L的离散信道；


x=generate_x(Nu+N0)+j*generate_x(Nu+N0);
y0=filter(ht,1,x)/norm(ht);

Yall=y0+sgma*(randn(size(y0))+j*randn(size(y0)));%加噪声
Y=Yall(N0:N0+Neu-1);
e=zeros(1,size(H,2));
e(1)=1;
w=e*H'*pinv(H*H');
wee=w/norm(w);
we=wee+randn(size(wee))*0.03;

%----------------------------恒模算法--------------------------------------------------
% w=zeros(1,N)*(1+j);
% w(1)=1+j;
w=we*(1+j);
mu=0.0002;
e_k=zeros(Neu,1);
for k=1:Neu-N0;
    x_k=Y(N0+k:-1:N0-N+1+k);
    y_k=w*x_k';
    %e_k=2*(y_k-y_k/abs(y_k));
    e_k(k)=(abs(y_k)).^2-2;
%w=w-mu*x_k*e_k';
    w=w-mu*e_k(k)*y_k*conj(x_k);
end;

S_e=zeros(1,Neu-N0);
for k=1:Neu-N0;
    x_k=Y(N0+k:-1:N0-N+1+k);
    S_e(k)=w*x_k';
    %e_k=2*(y_k-y_k/abs(y_k));
end;

figure(1),plot(Y,'*'),title('均衡前接收到的信号'),axis([-2.5 2.5 -2.5 2.5]);
figure(2),plot(S_e,'*'),
title('均衡后接收到的信号');


%--------------------------------------------------------------------------------------------------
function [X]=generate_x(n)
% use the matlab command "rand" to generate the binary,equally and independently distributed input data 
rX=rand(1,n);
tX=(rX>=0.5);
tX=2*tX-1;
X=tX;

%---------------------------------------------------------------------------------------
function [N_gau,seedn]=formNoise(Y,p)
L=size(Y,1)/p;
Ntemp=randn(p,size(Y,2));
N_gau=Ntemp;
Tr2=diag(ones(size(Y,2)-1,1),1);
Tr2(size(Y,2),1)=1;
for mm=1:L-1
    Ntemp=Ntemp*Tr2;
    Ntemp(:,1)=randn(p,1);
    N_gau=cat(1,N_gau,Ntemp);
end
N_gau=N_gau/norm(N_gau);