clear all;
N=5000;
M=31;


d=1;
M1=2;

   %sgma=sqrt()
   for j=1:N
       temp=rand;
       dsource(j)=1+floor(M1*temp);
   end;
mapping=[d*i
            
            -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,15,'measured');
figure(2);
plot(ChSym,'.b');
axis([-6.0 6.0 -6.0 6.0]);

% MMMA均衡
delt=0.0001;
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:1000
    
  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;
M1=[i -i];
ss=y;     
 figure(6);
  plot(ss,'b.');
  hold on
   plot(M1,'ro');
   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