% function Y_Y=BPSK_myself_VI()
%用自己編的卷積碼編解函數(shù)cnv_encd、viterbi仿真OFDM 
clear;


% tic
N=64; %The number of carriers

% h=gausComplex(1,3);
h=[1];
%h=[0.407,0.815,0.407];


L=length(h);
v=length(h)-1;%the size of the cyclic prefix


b=1;%For QPSK, one symbol carries  source bits



CBPSK=[-1,1];
BBPSK=[0;1];

   
srcBlockSize=N*b;
q=2^b;%the size of the signal set


Es=real(sum(CBPSK.*conj(CBPSK))/q);%ES= sum((abs(C16QAM)).^2)/qq;


Es=Es/N;% If noise is generated in time domain, Es should be multiplied by a factor 1/N



Eb=Es/b;

Na=N+v;%The size of the transmitted block with cyclic prefix

Nb=N+v+L-1;%The size of the liner convolutional result of the transmitted block and the channel

noiseVec=zeros(1,Nb);
X=zeros(1,N); % one source symbol block
x=zeros(1,N); % for the IFFT of X
xt=zeros(1,Na); % transmitted time domain block with cyclic prefix
yt=zeros(1,Nb);%received time domain block

y=zeros(1,N);%received time domain signal removal of CP
Y=zeros(1,N);%the vector fot FFT of y

XD=zeros(1,N);
dist=zeros(1,q);
src=zeros(1,srcBlockSize);
rev=zeros(1,srcBlockSize);
index=0;

H=fft(h,N); %the frequency domain character of the channel


randn('state',sum(100*clock));


recIndex=1;%an index for recoder buffer

for EbN0=0:1:0
    fprintf('%d>',EbN0)

    N0=Eb*10^(-EbN0/10);
    noiseRoot=sqrt(N0/2);
        
    errorCount=0;
    testLength=0;
    BER=0;
    while(1)
        testLength=testLength+srcBlockSize;
        input=(randn(1,srcBlockSize)>0);% Generate a source block (equal probability)
 
 %對(duì)原信息編碼 
k00=1;
G=[1 1 1;1 0 1];
src=cnv_encd(G,k00,input);
[row_src col_src]=size(src) ;    
              
        %The following part complement the conversion of binary to q-QAM symbol
        %a faster method is to use an index table,here we convert binary to symbol by calculating
        for k=1:1:N
            startIndex=(k-1)*b+1;
            endIndex=startIndex+b-1;
            vec=src(startIndex:endIndex);
            for t=1:1:q
                temp=sum(abs(vec-BBPSK(t,:)));
                if(temp<1e-5)
                    index=t;
                    break;
                end
            end
            X(k)=CBPSK(index);
        end
        
        %Mapping for frequency domain to time domain (multi-carrier modulation)
        x=ifft(X); %IFFT
        
        %Add CP
        xt(1:v)=x( (N-v+1):N);
        xt((v+1):(N+v))=x(1:N);
        
        %Passing the channel with noise
        yt=conv(xt,h);
        noiseVec=randn(1,Nb)+j*randn(1,Nb);%Generate noise
        noiseVec=noiseVec*noiseRoot;
        yt=yt+noiseVec;
        
        %Remove the cyclic prefix
        y=yt((v+1):(v+N));        
        
        %FFT of y
        Y=fft(y);
        
        %Equalizer on the frequency domain
        XD=Y./H;

        for k=1:1:N
            %minminum distance decision
            dist=abs(XD(k)-CBPSK);
             [temp,index]=min(dist);
             vec=BBPSK(index,:);
             startIndex=(k-1)*b+1;
             endIndex=startIndex+b-1;
             rev(startIndex:endIndex)=vec;
             
        end
        
        [row_rev col_rev]=size(rev);
        decode=viterbi(G,k00,rev);
        [row_input col_input]=size(input);
        [row_decode col_decode]=size(decode);
        
        
        errorNum=sum(xor(input,rev));
                
        errorCount=errorCount+errorNum;
        BER=errorCount/testLength;
        
        
                
        
        % plot(real(XD),imag(XD),'*'); %This can be used to watch the received signal set
        
        if(BER>0)
            if(EbN0<=9)
                temp=200/BER;
            else
                temp=100/BER;
            end
                       
            if(testLength>temp) 
                if(testLength>150000)
                    break;    
                end
            end
            
        end
        if(testLength>900000000)
%         if(testLength>90000000)
            break;    
        end
        
    end
    BER_Rec(recIndex)=BER;
    EbN0_Rec(recIndex)=EbN0;
    testLength_Rec(recIndex)=testLength;
    recIndex=recIndex+1;
end

%Save the result
% save BPSK_myself_myself_db0to8 BER_Rec;
BER_Rec
%Watch the results
semilogy(EbN0_Rec,BER_Rec,'k-s');

% time_in_hour=toc/3600
% save data\BPSK_myself_time time_in_hour