clear all; close all
Bit_num=2000;
 
 Fc=1800;
 Fs=18000;
 snr=5;
 Ts=0.001;
 t=0:1/Fs:Ts-1/Fs;
 LenSymbol=length(t);
 BW=2/Ts;
[b,a]=butter(5,2*BW/Fs);
Singal_emit=[];
Source=randint(1,Bit_num,[],0);%%信源

 
  TEMP=1-2*Source;%%%%%%%%%%%%%%調(diào)制
        for i_num=1:Bit_num
            Singal_emit=[Singal_emit TEMP(i_num)*cos(2*pi*Fc*t)];
        end
%         plot(Singal_emit(1:length(t)));
% m=0;
% for snr=-10:0.5:15
%     m=m+1;
        
        amp=sqrt(0.5/(10^(snr/10))); %%%%%%%%%加噪聲
Signal_rec1=Singal_emit+amp*randn(size(Singal_emit));

% plot(Signal_rec);
% vn=randn(1,T)+sqrt(-1)*randn(1,T);   % AWGN noise 
% vn=vn/norm(vn)*10^(-dB/20)*norm(x)

% plot(sigaddn);
Signal_rec=multipath(Signal_rec1);%%%加多徑
% Signal_rec=Singal_emit;

Base_wave_real=[];
Base_wave_imag=[];
 
           for i_num=1:Bit_num
            TEMP1=Signal_rec((i_num-1)*LenSymbol+1:i_num*LenSymbol);
            %             Xk=abs(fft(TEMP,4096));
            %             index=(0:round(2048-1));
            %             Xkf=Xk(index+1);
            %             indexf=find(Xkf==max(Xkf));
            %             Fc=indexf*Fs/4096;
            %             plot(index*Fs/4096,Xkf)
            %             pause
            %             periodogram(TEMP,[],4096,Fs);
            %               pause;
            TEMP_real=filter(b,a,TEMP1.*cos(2*pi*Fc*t));
%  TEMP_real=TEMP1.*cos(2*pi*Fc*t);
            TEMP_imag=filter(b,a,-1*TEMP1.*sin(2*pi*Fc*t));
            %              TEMP_real=filter(b,a,TEMP.*cn);
            %             TEMP_imag=filter(b,a,-1*TEMP.*sn);
            Base_wave_real=[Base_wave_real TEMP_real];
            Base_wave_imag=[Base_wave_imag TEMP_imag];
           end
        sample_real=[];
sample_imag=[];

        indexn=round(LenSymbol/2);
        for i_num=1:Bit_num
            TEMP_base_real=Base_wave_real((i_num-1)*LenSymbol+1:i_num*LenSymbol);
            TEMP_base_imag=Base_wave_imag((i_num-1)*LenSymbol+1:i_num*LenSymbol);
            sample_real_temp=TEMP_base_real(indexn);
            sample_imag_temp=TEMP_base_imag(indexn);
    %           sample_real_temp=sum(TEMP_base_real);
%               sample_imag_temp=sum(TEMP_base_imag);
            sample_real=[sample_real sample_real_temp];
            sample_imag=[sample_imag sample_imag_temp];
        end
        Sample=sample_real+j*sample_imag; 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求未均衡誤碼率
         g=0;
        for k=1:Bit_num
            if sample_real(k)<0
                newdata1(k)=1;
            else
                newdata1(k)=0;
            end
               if newdata1(k)==Source(k)
                   g=g+1;
               end
        end
        Pe1=(Bit_num-g)/Bit_num;
     
        %%%%%%%%%%均衡算法（LMS）
%         TEMP=TEMP+sqrt(-1)*TEMP;
        trseq=TEMP(1:800);
        trainlen=length(trseq);
        inseq=sample_real;
%       inseq=Sample;
        lev=40;
       wei=zeros(lev,1);
        miu=0.005;
        index=lev:-1:1;
        newinseq=[inseq,zeros(1,lev-1)];
        scatequ=[]; 
 for n=1:trainlen
    u=newinseq(index);
    d(n)=trseq(n);
    e(n)=d(n)-wei'*u.';
    mse(n)=abs(e(n))^2;
       wei=wei+miu*u.'*e(n)';
    index=index+1;
end
snrout=[];
 
for n=trainlen+1:length(inseq)
    u=newinseq(index);
    y=wei'*u.';
    scatequ=[scatequ,y];
      snrout=[snrout (abs(TEMP(n)-y))^2];
    [pan,pan_equ]=pan_2psk(y);
    d(n)=pan_equ;
       e(n)=d(n)-y;
    mse(n)=abs(e(n))^2;
      wei=wei+miu*u.'*e(n)';
    index=index+1;
end  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%均衡后判決      
        g=0;
        for k=1:Bit_num-trainlen
            if scatequ(k)<0
                newdata(k+trainlen)=1;
            else
                newdata(k+trainlen)=0;
            end
               if newdata(k+trainlen)==Source(k+trainlen)
                   g=g+1;
               end
        end
        Pe2=(Bit_num-g-trainlen)/(Bit_num-trainlen)
% end
c=trainlen+1:length(inseq);
d=trainlen+1:length(inseq);
    
  plot(c,scatequ,'.');
% plot(c,sample_real(d));
xlabel('碼元個數(shù)');
ylabel('幅度');
% x=-10:0.5:15;
% plot(x,Pe1);
% hold on
% Pe3=Pe2-0.2;
% for j=1:length(x)
%     if Pe3(j)<0
%         Pe3(j)=0
%     end
% end
% plot(x,Pe3,'--');
% xlabel('信噪比/dB');
% ylabel('誤碼率Pe');


%  end
% x=-10:0.5:15;
% % plot(x,Pe1);
% % hold on
% plot (x,Pe2);
% xlabel('SNR/dB');
% ylabel('誤碼率Pe');
% title('psk在高斯白噪聲與多徑情況下的誤碼率性能');
% scatterplot(sample);