%-----------------------------------------------------------------
%Simulation File
%OFDM System with 16QAM Modulation
%Self InterCarrier Interference Cancellation
%Frequency Offset Estimation Algorithm Implemented
%Size of FFT/IFFT=16
%Needed M-Files
%  QAMlev16
%  QAM16
%  gngauss
%  discrm16
%  revsymQAM16
%  levtobin16
%Author: 王博
%Filename: selficiawgnofdm.m
%-----------------------------------------------------------

%clearing and setting variables in the MATLAB environment
clear
home
format short g
N=16;                       %size of FFT
iterate=10000;                %number of iterations per Eb/No snr value
error=0:2:30;                 %Eb/No values to be scanned
%declaration of frequency offset corresponding from 0% to 10%
%with 2% increments 
%1% computed as (1/(size of FFT))/100
%so for FFT size 16 --> 1%=1/(16*100)=0.000625
freqerr=0:.00125:.00625;     %frequency offset with 2% increments
freqres=[];                 %matrix:BER values will be saved
echo off
%there will be 3 loops for the whole simulation
%for each frequency offset value the snr is sweeped for each snr value there is an repetitive iteration

for ff=1:length(freqerr) %first loop   (frequency offset)
for dif=1:length(error) %second loop   (snr values)

echo off;
count=0;
eroorcount=0;      %variable that monitors the number of bit
                 %errors

%computation of the corresponding sgma or variance of the
%Gaussian Noise generator. The whole system is normalized meaning
%the symbol (QAM16) signaling is kept constant while increasing the 
%variance of the GN Generator
d=1;
Eav=10*d^2;                
%average power of a 16 QAM system
snr=10^(error(diff)/10);    %computation of linear snr; error
                           %is in dB
sgma=sqrt(16*Eav/(8*snr)); %computation of the variance of
                       %Gaussian Generator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Numerous iteration starts here; third loop
for go=1:iterate

%start of transmitter model
%generate only half of the input bits of the conventional ofdm
y=rand(1,32);     %generates random values between 0 to 1
o=y>.5;          %floors or ceils; converts to logic 1 or 0

%also half the number of level conversions
for B=0:7          %groups the bits into sets of 4 numbers
sym(B+1, 1:4)=o(4*B+1:4*B+4);   %divides into 4 bit chunks
end

%QAM leveler
for B=0:7           %assigns a level number per combination of 4 numbers
symlev(B+1)=QAMlev16(sym(B+1, 1:4));
%QAMlev16 returns the level number
end

%generation of 16 QAM mapped symbols
for B=1:8          %each level will be assigned to a QAM constellation
aftQAM16(B)=QAM16(symlev(B));  %maps to a 16 point QAM point
end

%implenmentation of Linear Self ICI Cancellation in the
%transmitter side
aftQAM16B=[];
for B=1:8  %each QAM symbol is interleaved with its’ negated value
   aftQAM16B=[aftQAM16B aftQAM16(B) -1*(aftQAM16(B))];
   %resulting vector has an original vector length
end

%computation of the IFFT; complex modulation
OUTIFFT=N.*ifft(aftQAM16B,N);
%end of transmitter model

%channel model starts here
for B=1:16, [a b]=gngauss(sgma);   %generate Gaussian Noise
%here we insert noise in the channel
   OUTIFFTplusnoisel(B)=OUTIFFT(B)+a+i*b;
%frequency offset is also inserted into the signal
  %resulting vector OUT(B) is now affected by frequency offset and Gaussian Noise
  OUT(B)=OUTIFFTplusnoisel(B)*exp(i*2*pi(B-1)*freqerr(ff));
end;
%end of channel model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%start of receiver model
%computation of the FFT; complex demodulation
OUTFFT1=fft(OUT,N)./N;

%take 2 received signals and take fft and average before
%getting the symbol comparison
%Self ICI cancellation demodulation algorithm shown here
%sort out the odd and even frequency bins
k=1:2:16;    %odd pointers
k2=2:2:16;   %even pointers
OUTFFTa=OUTFFT1(k);   %odd frequency bin symbols
OUTFFTb=OUTFFT1(k2);  %even frequency bin symbols
%ICI Cancellation takes place right here
OUTFFTeven=(OUTFFTa-OUTFFTb0)/2;
%summation of two negated similar symbols and averaged will
%correspond to the symbol generated by the QAM16 constellation
%mapper only half of the symbols are seen in here

%we calculate the distances of the OUTFFT points to possible
%QAM constellation points
%the one with the least distances will assign as output
for B=1:8
   outsym(B)=discrim16(OUTFFTeven(B)); %return the nearest constellation coordinates
end

%reassigns the constellation coordinate to a voltage level
for B=1:8
revsymlev16(B)=revsymQAM16(outsym(B));
%voltage level corresponds to the assignment to QAMlev16
end

%reverse of QAM leveler
for B=1:8
outbin(B,1:4)=levtobin16(revsvmlev16(B));
%reassigns the voltage level to a binary pair
end
%end of receiver model

%start of bit error rate computation
%evaluation of the bit discrepancies found per 4 bit sets
for B=1:8
if outbin(B,1:4)==sym(B,1:4)
   errorcount=errorcount;     %perfect reception is achieved
else
   difference=sum(xor(outbin(B,1:4)));
%count the bit differences
errorcount=errorcount+difference;
%accumulate the bit errors
end          %end for the if statement
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end          %back to the numerous iteration
             %corresponds to the for loop of “for go=1:iterate”

%after counting all the bit errors per set of parameters
%compute the normalized number of errors accumulated
%with respect to the number of bits generated
bang(diff)=errorcount/(iterate*8*4);

end            %now update for new snr value
               %corresponds to for loop “for diff=1:length(error)”
               second loop (snr values)

%save the error count values for all levels of frequency offset
freqres=[freqres;bang];

end         %now update for a new frequency offset value
            %corresponds to for loop ff=1:length(freqerr)

a=clock;    %the time information
%saving the significant information for this simulation run
save AWGNselficiqwgnofdm freqres error N iterate d freqerr a 
%displaying the BER curves
%figure; semilogy(error,freqres);