mma_64.m

MMA盲均衡算法

clear all;
N=5000;
M=31;


d=1;
M1=64;

   %sgma=sqrt()
   for j=1:N
       temp=rand;
       dsource(j)=1+floor(M1*temp);
   end;
mapping =  [-7*d+7*d*i
            -5*d+7*d*i
            -3*d+7*d*i
            -1*d+7*d*i
             1*d+7*d*i
             3*d+7*d*i
             5*d+7*d*i
             7*d+7*d*i
            -7*d+5*d*i
            -5*d+5*d*i
            -3*d+5*d*i
            -1*d+5*d*i
             1*d+5*d*i
             3*d+5*d*i
             5*d+5*d*i
             7*d+5*d*i
            -7*d+3*d*i
            -5*d+3*d*i
            -3*d+3*d*i
            -1*d+3*d*i
             1*d+3*d*i
             3*d+3*d*i
             5*d+3*d*i
             7*d+3*d*i
            -7*d+1*d*i
            -5*d+1*d*i
            -3*d+1*d*i
            -1*d+1*d*i
             1*d+1*d*i
             3*d+1*d*i
             5*d+1*d*i
             7*d+1*d*i
            -7*d-1*d*i
            -5*d-1*d*i
            -3*d-1*d*i
            -1*d-1*d*i
             1*d-1*d*i
             3*d-1*d*i
             5*d-1*d*i
             7*d-1*d*i
            -7*d-3*d*i
            -5*d-3*d*i
            -3*d-3*d*i
            -1*d-3*d*i
             1*d-3*d*i
             3*d-3*d*i
             5*d-3*d*i
             7*d-3*d*i
            -7*d-5*d*i
            -5*d-5*d*i
            -3*d-5*d*i
            -1*d-5*d*i
             1*d-5*d*i
             3*d-5*d*i
             5*d-5*d*i
             7*d-5*d*i
            -7*d-7*d*i
            -5*d-7*d*i
            -3*d-7*d*i
            -1*d-7*d*i
             1*d-7*d*i
             3*d-7*d*i
             5*d-7*d*i
             7*d-7*d*i];
         for j=1:N
     SymQAM(j)=mapping(dsource(j),1);
 end;

figure(1);
plot(SymQAM,'.b');

Ch= [0.2 0.5 1 -0.1];
SymCh= filter(Ch,1,SymQAM);
ChSym = awgn(SymCh,25,'measured');
figure(2);
plot(ChSym,'.b');
axis([-30.0 30.0 -30.0 30.0]);

% MMMA均衡
delt=0.000001;
R=13.2;
%b=0.00000001;
%A_r=mean(abs(real(SymQAM)).^2);
%B_r=mean(abs(real(SymQAM)).^4);
%R_r=B_r/A_r;
%A_i=mean(abs(imag(SymQAM)).^2);
%B_i=mean(abs(imag(SymQAM)).^4);
%R_i=B_i/A_i;
%R=13.2;

est_c=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1+0j 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
est_cr=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.01 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
est_ci=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.01 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
C=7;

for k=1:200
    
  for j=M:N
      
    y(j-M+1)=ChSym(j:-1:j-M+1)*conj((est_c)');
    y_r(j-M+1)=real(y(j-M+1));
    y_i(j-M+1)=imag(y(j-M+1));
    est_cr=est_cr-delt*( y_r(j-M+1)*( y_r(j-M+1).^2-real(SymQAM(j-M+1)).^2))*conj(ChSym(j:-1:j-M+1));
    est_ci=est_ci-delt*( y_i(j-M+1)*( y_i(j-M+1).^2-imag(SymQAM(j-M+1)).^2))*conj(ChSym(j:-1:j-M+1));
    est_c=est_cr+i*est_ci;
   %e(j-M+1)=abs(y(j-M+1)).^2-R;
   % mse2(j-M+1)=e(j-M+1).^2;
   % est_c=est_c-delt*e(j-10)*y(j-10)*conj(ChSym(j:-1:j-10));
     %if(j>31)
     %   MSE2=mean(mse2(j-10:-1:j-30));
     %  MSE1=mean(mse2(j-11:-1:j-31));
        
     % MSE2(i)=mean(mse2);
     %delt=delt+b*(MSE2-MSE1);
     %MSE1=MSE2;
     %end;

  end;
end;
M6=[-7+7*i -5+7*i -3+7*i -1+7*i 1+7*i 3+7*i 5+7*i 7+7*i -7+5*i -5+5*i -3+5*i -1+5*i 1+5*i 3+5*i 5+5*i 7+5*i -7+3*i -5+3*i -3+3*i -1+3*i 1+3*i 3+3*i 5+3*i 7+3*i -7+i -5+i -3+i -1+i  1+i  3+i 5+i 7+i  -7-i -5-i -3-i -1-i 1-i 3-i 5-i 7-i -7-3*i -5-3*i -3-3*i -1-3*i 1-3*i  3-3*i  5-3*i 7-3*i -7-5*i -5-5*i -3-5*i  -1-5*i 1-5*i 3-5*i  5-5*i 7-5*i -7-7*i -5-7*i -3-7*i -1-7*i 1-7*i 3-7*i 5-7*i 7-7*i];
ss=y;     
 figure(6);
  plot(ss,'b.');
  hold on
   plot(M6,'ro');
    
  
  

%figure(4);
%sin=1:N/4-10; 
%plot(sin,mse2,'R');
%grid;

%山峰聚类
ss=y;     
 figure(6);
  plot(ss,'b.');
ss1=ss/max(abs(ss));
for j=1:6
   for k=1:6
    V(1,(j*6-6+k))=-1+(j-1)*0.4+i*(-1+(k-1)*0.4);
   end
end
a=0.0707;
b=0.1061;

 A=1/(2*(a.^2));
 B=1/(2*(b.^2));
% ra=0.2;
% rb=0.3;


% A=4/(ra.^2);
% B=4/(rb.^2);
n=length(ss1);
n2=length(V);

p=zeros(n2,1);

for i=1:n2

    for j=1:n
        p(i)=p(i)+exp(-A*((abs(V(i)-ss1(j))).^2));
    end
end
[pk(1,1),m]=max(p);
x(1,1)=V(m);

for i=1:n2
        p(i)=p(i)-pk(1,1)*exp(-B*((abs(V(i)-x(1,1))).^2));
    end

[pk(2,1),m]=max(p);
x(2,1)=V(m);
j=2;
while  pk(j,1)>=0.05*pk(1,1)
    
    
    for i=1:n2
        p(i)=p(i)-pk(j,1)*exp(-B*((abs(V(i)-x(j,1))).^2));
    end
    j=j+1;
    [pk(j,1),m]=max(p);
    x(j,1)=V(m);
end
j=j-1;
ss2=x(1:j,1);
   figure(3);
  plot(ss1,'b.');
  hold on
   plot(ss2,'ro');
axis([-1.0 1.0 -1.0 1.0]);



% 减法聚类
ss=y;     
 figure(6);
  plot(ss,'b.');
ss1=ss/max(abs(ss));
  
ra=0.2;
rb=0.3;
n=length(ss1);

A=4/(ra.^2);
B=4/(rb.^2);
p=zeros(n,1);

for i=1:n

    for j=1:n
        p(i)=p(i)+exp(-A*((abs(ss1(i)-ss1(j))).^2));
    end
end
[pk(1,1),m]=max(p);
x(1,1)=ss1(m);

    for i=1:n
        p(i)=p(i)-pk(1,1)*exp(-B*((abs(ss1(i)-x(1,1))).^2));
    end

[pk(2,1),m]=max(p);
x(2,1)=ss1(m);
j=2;
while  pk(j,1)>=0.15*pk(1,1)
    
    
    for i=1:n
        p(i)=p(i)-pk(j,1)*exp(-B*((abs(ss1(i)-x(j,1))).^2));
    end
    j=j+1;
    [pk(j,1),m]=max(p);
    x(j,1)=ss1(m);
end
j=j-1;
ss2=x(1:j,1);
   figure(3);
  plot(ss1,'b.');
  hold on
   plot(ss2,'ro');
% Ss1=abs(ss2)/max(abs(ss2));
%   
% ra=0.20;
% rb=0.30;
% N=length(Ss1);
% 
% A=4/(ra.^2);
% B=4/(rb.^2);
% P=zeros(N,1);
% 
% for i=1:N
% 
%     for j=1:N
%         P(i)=P(i)+exp(-A*((abs(Ss1(i)-Ss1(j))).^2));
%     end
% end
% [Pk(1,1),M]=max(P);
% X(1,1)=Ss1(M);
% Xx(1,1)=ss2(M);
%     for i=1:N
%         P(i)=P(i)-Pk(1,1)*exp(-B*((abs(Ss1(i)-X(1,1))).^2));
%     end
% 
% [Pk(2,1),M]=max(P);
% X(2,1)=Ss1(M);
% Xx(2,1)=ss2(M);
% j=2;
% while  Pk(j,1)>=0.15*Pk(1,1)
%     
%     
%     for i=1:N
%         P(i)=P(i)-Pk(j,1)*exp(-B*((abs(Ss1(i)-X(j,1))).^2));
%     end
%     j=j+1;
%     [Pk(j,1),M]=max(P);
%     X(j,1)=Ss1(M);
%     Xx(j,1)=ss2(M);
% end
% j=j-1;
% 
% Ss2=Xx(1:j,1);
% NN=length(Ss2);
% X=X(1:NN,1);
% [AA,bb]=max(X);
% j=1;
% for i=1:N
%     if min(abs(abs(Ss1(i))*ones(NN,1)-X))==abs(abs(Ss1(i))-X(bb));
%         ss3(j,1)=ss(i);
%         j=j+1;
%     end
% end
% NNN=j-1;
% j=1;
% for i=1:NNN
%     
% if angle(ss3(i))>=0
%     ss4(j,1)=ss3(i);
%     j=j+1;
% end
% end
% 
% qq=min(angle(ss4))-pi/4;
% qq1=min(angle(ss4));
% if NNN==4
%     ss2=ss2*exp(-qq*((-1).^0.5));
% else
%       ss2=ss2*exp(-qq1*((-1).^0.5));
%   end
%      figure(4);
%   plot(ss2,'b.');
%   
j=1;
for i=1:16
    if and(real(ss2(i))>=0,imag(ss2(i))>=0)
        ss3(j,1)=ss2(i)
        j=j+1;
    end
end
j=j-1;
[qq,iq]=max(abs(ss3));
qq1=angle(ss3(iq))-pi/4;
ss4=ss2*exp(-qq1*((-1).^0.5));
 figure(4);
 plot(ss4,'b.');
 
 
  xxx=48/(max(imag(ss4))-min(imag(ss4)));
    ss5=ss4*xxx;
    figure(5);
 plot(ss5,'b.');
 axis([-32 32 -32 32]);