%OFDM程序示例
%信道選用 J.G.Proakis Digital Communications 影印版 第四版 P631
%采用64子信道，應用16QAM傳輸，每個子載波均傳輸相同的比特數4。
%不使用信道編碼；使用迫零均衡和MMSE均衡，比較兩者性能。
%試應用仿真的方法畫出性能曲線。

close all;
clear all;


N=64; %The number of carriers


%h=[1];
%h=[0.407,0.815,0.407];
h=[0.04,-0.05,0.07,-0.21,-0.5,0.72,0.36,0,0.21,0.03,0.07];


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


b=4; %For 16QAM, one symbol carries  source bits

C16QAM=[3+3j,1+3j,-1+3j,-3+3j,-3+j,-1+j,1+j,3+j,3-j,1-j,-1-j,-3-j,-3-3j,-1-3j,1-3j,3-3j];
B16QAM=[0 0 0 0;0 0 0 1;0 0 1 1;0 0 1 0;0 1 1 0;0 1 1 1;0 1 0 1;0 1 0 0;
        1 1 0 0;1 1 0 1;1 1 1 1;1 1 1 0;1 0 1 0;1 0 1 1;1 0 0 1;1 0 0 0];
srcBlockSize=N*b;

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

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

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);
YD=zeros(1,N);
dist=zeros(1,q);
dist1=zeros(1,q);
src=zeros(1,srcBlockSize);
rev=zeros(1,srcBlockSize);
index=0;
index1=0;

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


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


recIndex=1;%an index for recoder buffer
recIndex1=1;
tic

for EbN0=0:1:20


    N0=Eb*10^(-EbN0/10);
    noiseRoot=sqrt(N0/2);
        
    errorCount=0;
    errorCount1=0;
    testLength=0;
    BER=0;
    BER1=0;
    while(1)
        testLength=testLength+srcBlockSize;
        src=(randn(1,srcBlockSize)>0);% Generate a source block (equal probability)
    
              
        %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-B16QAM(t,:)));
                if(temp<1e-5)
                    index=t;
                    break;
                end
            end
            X(k)=C16QAM(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;  %Zero forcing equalizing
        YD= Y.*((conj(H))./((H.*conj(H))+mean(noiseVec.*conj(noiseVec))/mean(x.*conj(x))));%MMSE forcing equalizing

        for k=1:1:N
            %minminum distance decision
            dist=abs(XD(k)-C16QAM);
            dist1=abs(YD(k)-C16QAM);
             [temp,index]=min(dist);
             [temp,index1]=min(dist1);
             vec=B16QAM(index,:);
             vec1=B16QAM(index1,:);
             startIndex=(k-1)*b+1;
             endIndex=startIndex+b-1;
             rev(startIndex:endIndex)=vec;
             rev1(startIndex:endIndex)=vec1;
             
        end
        errorNum=sum(xor(src,rev));
        errorNum1=sum(xor(src,rev1));
        errorCount=errorCount+errorNum;
        errorCount1=errorCount1+errorNum1;
        BER=errorCount/testLength;
        BER1=errorCount1/testLength;
        
                
        
        % plot(real(XD),imag(XD),'*'); %This can be used to watch the received signal set
        
        if((BER>0)&&(BER1>0))
            if(EbN0<=9)
                temp=200/BER;
                temp1=200/BER1;
            else
                temp=100/BER;
                temp1=100/BER1;
            end
                       
            if((testLength>temp)&&(testLength>temp1))
                if(testLength>150000)
                    break;    
                end
            end
            
        end
        
        if(testLength>90000000)
            break;    
        end
        
    end
    BER_Rec(recIndex)=BER
    BER_Rec1(recIndex)=BER1
    EbN0_Rec(recIndex)=EbN0
    testLength_Rec(recIndex)=testLength;
    recIndex=recIndex+1;
  
end

toc
%Save the result
save ofdm_data EbN0_Rec;
save ofdm_data BER_Rec;
save ofdm_data BER_Rec1;

%Watch the results
semilogy(EbN0_Rec,BER_Rec,'r-');
xlabel('Eb/N0 in dB');
ylabel('Bit error rate');
hold on
semilogy(EbN0_Rec,BER_Rec1,'*-');
xlabel('Eb/N0 in dB');
ylabel('hard decoding bit error rate');
hold on
legend('Zero forcing equalizing','MMSE equalizing');
grid