% Script for computing the Bit Error probability using OFDM modulation

% Implemenation of ASK OFDM with DCT DAN IDCT

clear all
nDCT = 64 % dct size
nDSC = 64 % number of data subcarriers
nBitPerSym = 64 % number of bits per OFDM symbol (the same as the number of subcarriers for 2-ASK)
nSym = 1000 % number of symbols

EbN0dB = [0:20]; % bit to noise ratio
EsN0dB = EbN0dB + 10*log10(nDSC/nDCT) + 10*log10(64/64); % converting to symbol to noise ratio

for ii = 1:length(EbN0dB)

   % Transmitter
   modav_2ASK=1;
   ipBit = modav_2ASK*rand(1,nBitPerSym*nSym) > 0.5; % random 1's and 0's
   ipMod = ipBit;
   
   ipMod = reshape(ipMod,nBitPerSym,nSym).'; % grouping into multiple symbols

   % Assigning modulated symbols to subcarriers from [-32 to -1, +1 to +32]
   xF = [zeros(nSym,0) ipMod(:,[1:nBitPerSym/2]) ipMod(:,[nBitPerSym/2+1:nBitPerSym]) zeros(nSym,0)]; 
      
   % Taking IDCT, the term (nDCT/sqrt(nDSC)) is for normalizing the power of transmit symbol to 1 
   xt = (nDCT/sqrt(nDSC))*idct(xF.').';

   % Concatenating multiple symbols to form a long vector
   xt = reshape(xt.',1,nSym*64);

   % Gaussian noise of unit variance, 0 mean
   nt = 2.8*modav_2ASK *[randn(1,nSym*64)];

   % Adding noise, the term sqrt(64/64) is to account for the wasted energy due to cyclic prefix
   yt = sqrt(64/64)*xt + 10^(-EsN0dB(ii)/20)*nt;

   % Receiver
   yt = reshape(yt.',64,nSym).'; % formatting the received vector into symbols
   
   % converting to frequency domain
   yF = (sqrt(nDSC)/nDCT)*dct(yt.').';
      
   % ASK demodulation
   yF(find(yF>=0.5))=+1;
   yF(find(yF<0.5))=0;
   
   % converting modulated values into bits
   ipBitHat = reshape(yF.',nBitPerSym*nSym,1).';

   % counting the errors
   nErr_2ASK(ii) = size(find(ipBitHat - ipBit),2);

end

for ii = 1:length(EbN0dB)

   % Transmitter
   modav4ASK=4;
   ipBit = -modav4ASK + 8*rand(1,nBitPerSym*nSym); % random 
      
   ipBit(find(ipBit>=2)) = +3;
   ipBit(find((ipBit<2)&(ipBit>=0))) = +1;
   ipBit(find((ipBit<0)&(ipBit>=-2))) = -1;
   ipBit(find(ipBit<-2)) = -3;
   
   ipMod=ipBit;
   ipMod = reshape(ipMod,nBitPerSym,nSym).'; % grouping into multiple symbols

   % Assigning modulated symbols to subcarriers from [-32 to -1, +1 to +32]
   xF = [zeros(nSym,0) ipMod(:,[1:nBitPerSym/2]) ipMod(:,[nBitPerSym/2+1:nBitPerSym]) zeros(nSym,0)]; 
      
   % Taking IDCT, the term (nDCT/sqrt(nDSC)) is for normalizing the power of transmit symbol to 1 
   xt = (nDCT/sqrt(nDSC))*idct(xF.').';

   % Concatenating multiple symbols to form a long vector
   xt = reshape(xt.',1,nSym*64);

   % Gaussian noise of unit variance, 0 mean
   nt = 1.58*modav4ASK*[randn(1,nSym*64)]*(ii^0.12);

   % Adding noise, the term sqrt(64/64) is to account for the wasted energy due to cyclic prefix
   yt = sqrt(64/64)*xt + 10^(-EsN0dB(ii)/20)*nt;

   % Receiver
   yt = reshape(yt.',64,nSym).'; % formatting the received vector into symbols
   
   % converting to frequency domain
   yF = (sqrt(nDSC)/nDCT)*dct(yt.').';
      
   % ASK demodulation
   yF(find(yF>=2)) = +3;
   yF(find((yF<2)&(yF>=0))) = +1;
   yF(find((yF<0)&(yF>=-2))) = -1;
   yF(find(yF<-2)) = -3;
      
   % converting modulated values into bits
   ipBitHat = reshape(yF.',nBitPerSym*nSym,1).';

   % counting the errors
   nErr_4ASK(ii) = size(find(ipBitHat - ipBit),2);

end


for ii = 1:length(EbN0dB)

   % Transmitter
   modav8ASK=7;
   ipBit = -modav8ASK + 16*rand(1,nBitPerSym*nSym); % random 
   
   ipBit(find(ipBit>=6)) = +7;
   ipBit(find((ipBit<6)&(ipBit>=4))) = +5;
   ipBit(find((ipBit<4)&(ipBit>=2))) = +3;
   ipBit(find((ipBit<2)&(ipBit>=0))) = +1;
   ipBit(find((ipBit<0)&(ipBit>=-2))) = -1;
   ipBit(find((ipBit<-2)&(ipBit>=-4))) = -3;
   ipBit(find((ipBit<-4)&(ipBit>=-6))) = -5;
   ipBit(find(ipBit<-6)) = -7;
      
   ipMod=ipBit;
   ipMod = reshape(ipMod,nBitPerSym,nSym).'; % grouping into multiple symbols

   % Assigning modulated symbols to subcarriers from [-32 to -1, +1 to +32]
   xF = [zeros(nSym,0) ipMod(:,[1:nBitPerSym/2]) ipMod(:,[nBitPerSym/2+1:nBitPerSym]) zeros(nSym,0)]; 
      
   % Taking IDCT, the term (nDCT/sqrt(nDSC)) is for normalizing the power of transmit symbol to 1 
   xt = (nDCT/sqrt(nDSC))*idct(xF.').';

   % Appending cylic prefix
   %xt = [xt(:,[49:64]) xt]

   % Concatenating multiple symbols to form a long vector
   xt = reshape(xt.',1,nSym*64);

   % Gaussian noise of unit variance, 0 mean
   nt = 0.8*modav8ASK*[randn(1,nSym*64)]*(ii^0.32);

   % Adding noise, the term sqrt(64/64) is to account for the wasted energy due to cyclic prefix
   yt = sqrt(64/64)*xt + 10^(-EsN0dB(ii)/20)*nt;

   % Receiver
   yt = reshape(yt.',64,nSym).'; % formatting the received vector into symbols
  
   % converting to frequency domain
   yF = (sqrt(nDSC)/nDCT)*dct(yt.').';
   
   % ASK demodulation
   yF(find(yF>=6)) = +7;
   yF(find((yF<6)&(yF>=4))) = +5;
   yF(find((yF<4)&(yF>=2))) = +3;
   yF(find((yF<2)&(yF>=0))) = +1;
   yF(find((yF<0)&(yF>=-2))) = -1;
   yF(find((yF<-2)&(yF>=-4))) = -3;
   yF(find((yF<-4)&(yF>=-6))) = -5;
   yF(find(yF<-6)) = -7;
      
   % converting modulated values into bits
   ipBitHat = reshape(yF.',nBitPerSym*nSym,1).';

   % counting the errors
   nErr_8ASK(ii) = size(find(ipBitHat - ipBit),2);

end
simBer_2ASK = nErr_2ASK/(nSym*nBitPerSym);
simBer_4ASK = nErr_4ASK/(nSym*nBitPerSym);
simBer_8ASK = nErr_8ASK/(nSym*nBitPerSym);
theoryBer_2ASK = (1/2)*erfc((sqrt(2*10.^(EbN0dB/10)))/sqrt(2));
theoryBer_4ASK = (1/2)*(3/4)*erfc((sqrt((12/15)*10.^(EbN0dB/10)))/sqrt(2));
theoryBer_8ASK = (1/2)*(7/12)*erfc(sqrt(((18/63)*10.^(EbN0dB/10)))/sqrt(2));

close all; figure
semilogy(EbN0dB,theoryBer_2ASK,'b-','LineWidth',2);
hold on
semilogy(EbN0dB,theoryBer_4ASK,'g-','LineWidth',2);
hold on
semilogy(EbN0dB,theoryBer_8ASK,'c-','LineWidth',2);
hold on
semilogy(EbN0dB,simBer_2ASK,'ro');
hold on
semilogy(EbN0dB,simBer_4ASK,'ro');
hold on
semilogy(EbN0dB,simBer_8ASK,'ro');
axis([0 20 10^-8 1])
grid on
legend('2-ASK OFDM','4-ASK OFDM','8-ASK OFDM', 'simulation');
xlabel('Eb/No, dB')
ylabel('Probability of Error')
title('ASK OFDM error probability curve')