mma_16.m

MMA盲均衡算法

clear all;
N=20000;
M=31;
T=100;
BitsTx = floor(rand(1,N)*2);
SymQAMtmp = reshape(BitsTx,4,N/4).';
SymQAMtmptmp= bi2de(SymQAMtmp,'left-msb');
QAMTable = [-3-3*i -3-i -3+3*i -3+i -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];
SymQAM = QAMTable(SymQAMtmptmp+1);
figure(1);
plot(SymQAM,'.b');
axis([-6.0 6.0 -6.0 6.0]);
Ch= [0.2 0.5 1 -0.1];
SymCh= filter(Ch,1,SymQAM);
ChSym = awgn(SymCh,15,'measured');
figure(2);
plot(ChSym,'.b');
axis([-6.0 6.0 -6.0 6.0]);

% MMMA均衡
delt=0.00001;
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 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.1 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.1 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/4
      
    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;
M5=2*[-3+5*i -1+5*i 1+5*i 3+5*i -5+3*i -3+3*i -1+3*i 1+3*i 3+3*i 5+3*i -5+i -3+i -1+i 1+i 3+i 5+i -5-i -3-i -1-i 1-i 3-i 5-i -5-3*i -3-3*i -1-3*i 1-3*i 3-3*i 5-3*i -3-5*i -1-5*i 1-5*i 3-5*i];
ss=y;     
 figure(6);
  plot(ss,'b.');
   hold on
   plot(M5,'ro');
  
    
  
  

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

%山峰聚类

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]);

%识别
M1=[i -i];
M2=[1 i -1 -i];
M3=[cos(2*pi/8)+sin(2*pi/8)*i -cos(2*pi/8)+sin(2*pi/8)*i -cos(2*pi/8)-sin(2*pi/8)*i cos(2*pi/8)-sin(2*pi/8)*i 1 i -1 -i];
M4=[-3+3*i -1+3*i 1+3*i 3+3*i -3+i -1+i 1+i 3+i -3-i -1-i 1-i 3-i -3-3*i -1-3*i 1-3*i 3-3*i ];
M5=[-3+5*i -1+5*i 1+5*i 3+5*i -5+3*i -3+3*i -1+3*i 1+3*i 3+3*i 5+3*i -5+i -3+i -1+i 1+i 3+i 5+i -5-i -3-i -1-i 1-i 3-i 5-i -5-3*i -3-3*i -1-3*i 1-3*i 3-3*i 5-3*i -3-5*i -1-5*i 1-5*i 3-5*i];
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];
n=length(ss);
n1=length(M1);
r=0.7;
for k=1:n1

    for j=1:n
        p1(k,j)=abs(M1(k)-ss(j));
    end
end

m1=numel(p1);
p1=reshape(p1,1,m1);

i1=0;
for k=1:m1
  if (p1(k)<r)
     i1=i1+1;
  end
end

n2=length(M2);
for k=1:n2

    for j=k:n
        p2(k,j)=abs(M2(k)-ss(j));
    end
end
m2=numel(p2);
p2=reshape(p2,1,m2);

i2=0;
for k=1:m2
  if (p2(k)<(r/1.414))
     i2=i2+1;
  end
end

n3=length(M3);
for k=1:n3

    for j=1:n
        p3(k,j)=abs(M3(k)-ss(j));
    end
end

m3=numel(p3);
p3=reshape(p3,1,m3);

i3=0;
for k=1:m3
  if (p3(k)<(r/2))
     i3=i3+1;
  end
end

n4=length(M4);
for k=1:n4

    for j=1:n
        p4(k,j)=abs(M4(k)-ss(j));
    end
end

m4=numel(p4);
p4=reshape(p4,1,m4);

i4=0;
for k=1:m4
  if (p4(k)<(r/2))
     i4=i4+1;
  end
end

n5=length(M5);
for k=1:n5

    for j=1:n
        p5(k,j)=abs(M5(k)-ss(j));
    end
end

m5=numel(p5);
p5=reshape(p5,1,m5);
i5=0;
for k=1:m5
  if (p5(k)<(r/4))
     i5=i5+1;
  end
end

n6=length(M6);
for k=1:n6

    for j=1:n
        p6(k,j)=abs(M6(k)-ss(j));
    end
end

m6=numel(p6);
p6=reshape(p6,1,m6);
i6=0;
for k=1:m6
  if (p6(k)<(r/5.656))
     i6=i6+1;
  end
end


% 减法聚类
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]);