m=150;
pi=3.14159;
n=10;           
id=sqrt(-1);
numerr=0;                                                               %誤碼率
numerr1=0;
iteration_num=1;
for a=1:iteration_num
%-------------------------------------
%碼元產生
ran=randn(1,m);                  %隨機產生Bpsk的信號比特1或-1
data_sequence=zeros(1,m);
for i=1:m
  if ran(1,i)>0
      data_sequence(1,i)=1;
  else
      data_sequence(1,i)=0;
  end   
end;

%--------------------------------------------
%23訓練序列+數據
sequence=[+1 +1 -1 +1 -1 +1 +1 +1 +1 -1 +1 +1 -1 -1 +1 +1 +1 -1 +1 -1 +1 +1 +1];
safeguard1=[-1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
safeguard2=[-1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];

%差分編碼
M=4;          
mapping=[-3 3 -1 1];
data_sequence_de=coding(data_sequence,M,mapping);
Ik=data_sequence_de(:,1)';
Qk=data_sequence_de(:,2)';
data=Ik+Qk*id;

%sequence=[safeguard1 sequence data safeguard2];
sequence=[sequence data];
%--------------------------------------------

%-------------------------------------------------------
%波形形成,內插和調制
[I_wave,Q_wave]=wavemod(sequence);
I_wave=interp(I_wave,4);
Q_wave=interp(Q_wave,4);
fs=5*128000*4;               
fc=455000;
t=0:1/fs:(length(I_wave)-1)*1/fs;
mod_train=I_wave.*cos(2*pi*fc*t)-Q_wave.*sin(2*pi*fc*t);
%-------------------------------------------------------


%------------------------------------------------------
%經信道加噪聲和多徑
snr=20;
channel=[1 0  0.5*(0.996+0.087*id)];
channel=[ 1+0j 0.3-0.06*id 0.2+0.1*id];
channel=[1 0.2+0.1*id];
channel=[0.1 1 0.2];
[channel_out,N0]=multi_channel(mod_train,channel,snr);
channel_out=channel_out(21:length(mod_train)+20);
%channel_out=channel_out(1:length(mod_train));
%------------------------------------------------------
%解調和抽取
demod_out=demod1(channel_out);

demod_out_extract=zeros(1,length(demod_out)/4);
for i=1:length(demod_out)/4
    demod_out_extract(i)=demod_out(4*(i-1)+1);
end;

%------------------------------------------------------
%信道估值
refe_sequence=[ -1 +1 +1 +1 +1 -1 +1 +1 -1 -1 +1 +1 +1 -1 +1];          %本地參考序列


corr=zeros(23,5);
test=0;
for k=1:23
  for i=1:5
    select=zeros(1,15);  
    for j=1:15
      select(j)=demod_out_extract(5*(j-1+k-1)+i);
    end;
     corr(k,i)=(sum(refe_sequence.*(select)));
     if abs(corr(k,i))>test                                            %記錄最大的相關值（位同步點位置p,q ）
       test=abs(corr(k,i));
       p=k;
       q=i;
     end;
     if k==5&i==3
         de=select;
     end;
  end;
end;
 
G=[15 -1 -1 -1 3
   -1 15 -1 -1 3
   -1 -1 15 -1 3
   -1 -1 -1 15 3
   -1 -1 -1 3 15];
A=[corr(p-2,q) corr(p-1,q) corr(p,q) corr(p+1,q) corr(p+2,q)].';
W=inv(G)*A;


Tw=0.1;                   %門限
for i=1:length(W)         %大于門限的信道參數保留,小于門限的信道參數置0
    if abs(W(i))<Tw
        W(i)=0;
    end;
end;
 

refe_sequence1=[1 refe_sequence(1:14)];
refe_sequence2=[-1 1 refe_sequence(1:13)];
refe_sequence3=[1 -1 1 refe_sequence(1:12)];
refe_sequence4=[+1 +1 -1 +1 refe_sequence(1:11)];
error=sum(abs(refe_sequence*W(1)+refe_sequence1*W(2)+refe_sequence2*W(3)+refe_sequence3*W(4)+refe_sequence4*W(5)-de))/15; %重構誤差

%------------------------------------------------------


%------------------------------------------------------
%碼元中點取樣,中點在第q個樣點                                                
y_in=zeros(1,length(demod_out_extract)/5);
for i=1:length(demod_out_extract)/5
    y_in(i)=demod_out_extract(5*(i-1)+q);
end;
%------------------------------------------------------


%------------------------------------------------------
%均衡


h_estimated=[W(1) W(2) W(3) W(4) W(5)];
%h_estimated=[W(3) W(4) W(5)];
        
L=5;
 %-------------------------------------------------
 % 由信道參數計算均衡器系數
    
    g=zeros(L,L);
    for i=1:L,
      for j=i:L,
        for n=1:L-(j-i), 
           g(i,j)=g(i,j)+h_estimated(n)*conj(h_estimated(n+j-i));           %g(i,j)=h(0)*conj(h(0+j-i))+h(1)*conj(h(1+j-i))+...+h(L)*conj(h(L+j-i))
        end;
        if i==j
          g(i,j)=g(i,j);%+N0;
        end;
      end;
    end;

    for i=1:L,
      for j=1:i,  
        g(i,j)=g(j,i);
      end;  
    end;
    
    f=zeros(1,L);
    for i=1:L,
       f(i)=conj(h_estimated(L+1-i));                                       % f=conj[h(K1),.....,h(1),h(0)]
    end;
    c_estimated1=inv(g)*f.';                                                %計算均衡器前饋級系數
    
    %-------------------------------------------------------------
    % Gauss-sidel迭代求前饋級系數(DSP實現時用此法代替矩陣求逆)
    c1=0;
    c2=0;
    c3=0;
    c4=0;
    c5=0;
    for i=2:5
        c1=c1+1/g(1,1)*(f(1)-g(1,1)*c1-g(1,2)*c2-g(1,3)*c3-g(1,4)*c4-g(1,5)*c5);
        c2=c2+1/g(2,2)*(f(2)-g(2,1)*c1-g(2,2)*c2-g(2,3)*c3-g(2,4)*c4-g(2,5)*c5);
        c3=c3+1/g(3,3)*(f(3)-g(3,1)*c1-g(3,2)*c2-g(3,3)*c3-g(3,4)*c4-g(3,5)*c5);   
        c4=c4+1/g(4,4)*(f(4)-g(4,1)*c1-g(4,2)*c2-g(4,3)*c3-g(4,4)*c4-g(4,5)*c5);
        c5=c5+1/g(5,5)*(f(5)-g(5,1)*c1-g(5,2)*c2-g(5,3)*c3-g(5,4)*c4-g(5,5)*c5);
    end;
    %-------------------------------------------------------------
    
    K2=L;                                                                  %計算均衡器反饋級系數    
    c_estimated2=zeros(1,K2);
    for i=1:K2
        for j=i:K2-1
          c_estimated2(i)=c_estimated2(i)-c_estimated1(j+1)*conj(f(j-i+1)); %c(k)=-(c(-K1)*f(k+K1)+c(-K1+1)*f(k+K1-1)+....+c(0)*f(k))
      end;
    end;    
    %---------------------------------------------------------------
    
    %---------------------------------------------------------------
    %均衡
   
   
   iout=zeros(1,length(y_in)-L);
   x_test=zeros(1,length(y_in)-L);
   y_test=zeros(1,length(y_in)-L);
   x_test1=zeros(1,length(y_in)-L);
   y_test1=zeros(1,length(y_in)-L);
   se_test=zeros(1,length(y_in)-L);
   eq_out=zeros(1,length(y_in)-L);
   uneq_out=zeros(1,length(y_in)-L);
   eq_out(21)=1;
   uneq_out(23)=1;
   y_in(22)=1;
   y_in(23)=1;
   for b=22:length(y_in)-L
    
    for j=1:L
           yy_k(j)=y_in(b+L-j);                                             %前饋級數據  [y(b+K1),y(b+K1-1),...,y(b)]       
    end;
    for j=1:K2                             
         y_k(j)=eq_out(b-j);                                                 %反饋級數據  [I~(b-1),I~(b-2),...,I~(b-K2)]
     end;
    yy1=yy_k(1:L).';
    yy2=y_k(1:K2).';
    iout1(1,b)=((c_estimated1).')*yy1; 
    iout2(1,b)=((c_estimated2))*yy2;
   
    iout(1,b)=iout1(1,b)+iout2(1,b);                                         % I^(b) 均衡器的輸出
    if real(iout(1,b))>=0.8536
        eq_out(1,b)=1;
    elseif real(iout(1,b))>=0.3536
        eq_out(1,b)=0.7071;
    elseif real(iout(1,b))>=-0.3536
        eq_out(1,b)=0;
    elseif real(iout(1,b))>=-0.8536
        eq_out(1,b)=-0.7071;
    else
        eq_out(1,b)=-1;
    end;   
    if imag(iout(1,b))>=0.8536
        eq_out(1,b)=eq_out(1,b)+1*id;
    elseif imag(iout(1,b))>=0.3536
        eq_out(1,b)=eq_out(1,b)+0.7071*id;
    elseif imag(iout(1,b))>=-0.3536
        eq_out(1,b)=eq_out(1,b)+0*id;
    elseif imag(iout(1,b))>=-0.8536
        eq_out(1,b)=eq_out(1,b)+(-0.7071)*id;
    else
        eq_out(1,b)=eq_out(1,b)+(-1)*id;
    end;   
    
    e(b)=real(eq_out(1,b))*real(eq_out(1,b-1))+imag(eq_out(1,b))*imag(eq_out(1,b-1));
    f(b)=imag(eq_out(1,b))*real(eq_out(1,b-1))-real(eq_out(1,b))*imag(eq_out(1,b-1));
   
    if e(b)>0                                                        %字符檢測  I~(b)
        x_test(b)=1;
    elseif e(b)<0
        x_test(b)=0;
    end;
    if f(b)>0                                                        %字符檢測  I~(b)
        y_test(b)=1;
    elseif e(b)<0
        y_test(b)=0;
    end;
        

   if abs(data_sequence(2*(b-22)+1)-x_test(b))>0.1
       numerr=numerr+1;
   end;
   
   if abs(data_sequence(2*(b-22)+2)-y_test(b))>0.1
       numerr=numerr+1;
   end;
   end;
   
   
   for b=24:length(y_in)-L
   if real(y_in(1,b))>=0.8536
        uneq_out(1,b)=1;
    elseif real(y_in(1,b))>=0.3536
        uneq_out(1,b)=0.7071;
    elseif real(y_in(1,b))>=-0.3536
        uneq_out(1,b)=0;
    elseif real(y_in(1,b))>=-0.8536
        uneq_out(1,b)=-0.7071;
    else
        uneq_out(1,b)=-1;
    end;   
    if imag(y_in(1,b))>=0.8536
        uneq_out(1,b)=uneq_out(1,b)+1*id;
    elseif imag(y_in(1,b))>=0.3536
        uneq_out(1,b)=uneq_out(1,b)+0.7071*id;
    elseif imag(y_in(1,b))>=-0.3536
        uneq_out(1,b)=uneq_out(1,b)+0*id;
    elseif imag(y_in(1,b))>=-0.8536
        uneq_out(1,b)=uneq_out(1,b)+(-0.7071)*id;
    else
        uneq_out(1,b)=uneq_out(1,b)+(-1)*id;
    end;   
    
    
    e1(b)=real(uneq_out(1,b))*real(uneq_out(1,b-1))+imag(uneq_out(1,b))*imag(uneq_out(1,b-1));
    f1(b)=imag(uneq_out(1,b))*real(uneq_out(1,b-1))-real(uneq_out(1,b))*imag(uneq_out(1,b-1));
    if e1(b)>0                                                        %字符檢測  I~(b)
        x_test1(b)=1;
    elseif e(b)<0
        x_test1(b)=0;
    end;
    if f1(b)>0                                                        %字符檢測  I~(b)
        y_test1(b)=1;
    elseif e(b)<0
        y_test1(b)=0;
    end;
    
   if abs(data_sequence(2*(b-24)+1)-x_test1(b))>0.1
       numerr1=numerr1+1;
   end;
   
   if abs(data_sequence(2*(b-24)+2)-y_test1(b))>0.1
       numerr1=numerr1+1;
   end;
end
end

numerr=numerr/iteration_num;
numerr1=numerr1/iteration_num;
unequalized=y_in(24:73);
equalized=iout(24:70);
figure(1);
plot(real(unequalized),imag(unequalized),'o')
axis('square')
axis([-1.5 1.5 -1.5 1.5]);
figure(2);
plot(real(equalized),imag(equalized),'o')
axis('square')
axis([-1.5 1.5 -1.5 1.5]);
a=1;

w=0.995;
p2=eye(5)*0.004;
c2=[0 0 0 0 0];
e2=zeros(1,length(y_in)-L);
iout2=zeros(1,length(y_in)-L);
ks2=zeros(5,length(y_in)-L);
xx2=0;
yy2=zeros(5,1);
 for b=22:length(y_in)-L
    
    for j=1:L
           yy_k(j)=y_in(b+L-j);                                             %前饋級數據  [y(b+K1),y(b+K1-1),...,y(b)]       
    end;
   iout2(1,b)=conj(c2)*yy1; 
   e2(b)=data(b-21)-iout2(b);
   xx2=w+(conj(yy1).')*p2*yy1;
   yy2=p2*yy1;
   ks2(:,b)=yy2/xx2;
   p2=(p2-ks2(:,b)*(conj(yy1).')*p2)/w;
   c2=c2+ks2(:,b).'*conj(e2(b));
   mse2(i)=abs(e2(b))^2;
end
a=1;