estcorrawgnofdm.m

正交频分复用系统(OFDM)的频率偏移误差抑制技术仿真程序

%------------------------------------------------------------------
%Simulation File
%OFDM System with 16 QAM Modulation
%Frequency Offset Estimation and SymbolCorrection Algorithm
%Implemented
%Size of FFT/IFFT=16
%Needed M-Files
%  QAMlev16
%  QAM16
%  gngauss
%  discrm16
%  revsymQAM16
%  levtobin16
%Author:  王博
%Filename: estcorrawgnofdm.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 repetitve iteration

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

   echo off;
   count=0;
   errorcount=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

Y=rand(1,64);    %generates random values between 0 to 1
o=Y>.5;         %floors or ceils; converts to logic 1 or 0

for B=0:15       %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:15
%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:16
%each level will be assigned to the QAM constellation
   aftQAM16(B)=QAM16(symlev(B));
%maps to a 16 point QAM point
   end
%computation of the IFFT; complex modulation
OUTIFFT=N.*ifft(aftQAM16,N);

%retransmission of the IFFT values as dictated by
%the frequency offset estimation and correction algorithm
doubIFFT=[OUTIFFT OUTIFFT];
%note :  different from conventional OFDM
%end of transmitter model

%channel model starts here
for B=1:length(doubIFFT),
   [a,b]=gngauss(sgma);       %generate of Gaussian noise
%here we insert noise in the channel
OUTIFFTplusnoisel(B)=doubIFFT(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*pi*(B-1)*freqerr(ff));
end;
%end of channel model

%start of receiver Model with frequency offset correction technique
%start of Frequency Offset Estimation and Correction Block
%separate received symbols into two sets
  OUT1=OUT(1:N);          %first set
  OUT2=OUT(N+1:2*N);      %second set
%each the second set values are aligned to the original values
%the only difference is that in the second set frequency offset
%is present
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%computation of the frequency offset estimates coming from each
%aligned pair
aveest=1/(2*pi*N)*atan(sum(imag(OUT2.*conj(OUT1)))./sum(real(OUT2.*conj(OUT1))));

%generate a time vector
k=0:(2*N-1);
%generate a frequency offset complex modulation with a negated frequency offset estimate
revOUT=OUT.*exp(i*2*pi*k*(-aveest));    %correction vector
%symbol correction is done by multiplying the distorted symbols
%to the correction vector

%compute the FFT; complex demodulation
decOUT=fft(revOUT,(2*N))./(2*N);
%notice that the FFT here is twice the size
%so TX is N IFFT
%but RX is 2*N FFT

%get the even bin result values from the FFT
kk=1:2:32;
evalOUT=decOUT(kk);

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

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

%reverse of QAM leveler
for B=1:16
outbin(B, 1:4)=levtobin16(revsymlev16(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:16
if   outbin(B,1:4)== sym(B,1:4)
errorcount=errorcount;
%perfect reception is achieved
else
  difference=sum(xor(outbin(B,1:4),sym(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*16*4);

end          %now update for new snr value
             %corresponds to for loop “for diff=1:lengthe(error)”
             %second loop (snt 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 “for ff=1:length(freqerr)”

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