ofdm.m

在matlab中的ofdm的编程

%ofdm simulation
%-------------------------preparation part---------------------------------
para=52;           %parallel channel number
fftlen=64;         %FFT length
noc=53;            %number of carrier 
nd=6;              %number of one loop 's ofdm symbol

knd=1;             

ml=2;              %modulation level
sr=250000;         %symbol rate
br=sr*ml;          %bit rate of ber carrier
gilen=16;          %length of guard interval    

%------------------------fading initialization-----------------------------
tstp=1/sr/(fftlen+gilen);                       %time resolution
itau=[6];                                       %the number of symbols to be delayed
n0=[6];                                         %number of wave to generate fading  
itnd1=[1000];                                   %fading counter(the number of fading counter to skip)
now1=6;                                         %number of direct wave impluse the multiple wave)
fd=150;                                         %the biggest doppler frequency
flat=0;                                         %flat frequency or not 

%----------------------main loop part--------------------------------------
nloop=1000;
noe=0;
nod=0;
eop=0;
nop=0;

ich0=zeros(fftlen, nd);
qch0=zeros(fftlen, nd);
paradata=zeros(fftlen, nd);
middlevalue1=zeros(fftlen, nd);
middlevalue2=zeros(fftlen, nd);
%计算误码率-----------------------------------------------------------------
numberofpacket=1;
SNR_start=0;
SNR_end=25;
p=1;
BER=zeros(numberofpacket, SNR_end-SNR_start+1);
for SNR=SNR_start:SNR_end
    for count=1:numberofpacket
         for iii=1:nloop
   
             seldata=rand(1,para*nd*ml)>0.5;                                %data the serial digital
    
             paradata=reshape(seldata, para, nd*ml);                        %change the serial digital to parallel data 
    %QPSK 调制
             paradata2=paradata.*2-1;
              for x=1:nd
                  ich0((1:para),x)=paradata2((1:para),x*2-1);
                  qch0((1:para),x)=paradata2((1:para),x*2);
              end
    %CE(channel estimate) date generation----------------------------------
            kndata=zeros(1,fftlen);
            kndata0=2.*(rand(1,52)<0.5)-1;
            kndata(2:27)=kndata0(1:26);
            kndata(39:64)=kndata0(27:52);
            ceich=kndata;
            ceqch=zeros(1,64);
    %------------------------data mapping----------------------------------
    
           ich1=zeros(fftlen, nd);
           qch1=zeros(fftlen, nd);   
           ich1(2:27,:)=ich0(1:26,:);
           qch1(2:27,:)=qch0(1:26,:);
           ich1(39:64, :)=ich0(27:52,:);
           qch1(39:64, :)=qch0(27:52,:);
   
           ich2=[ceich.' ich1];
           qch2=[ceqch.' qch1];
         
    %IFFT------------------------------------------------------------------
           ich3=zeros(fftlen, nd+1);
           qch3=zeros(fftlen, nd+1);
           middlevalue1=ich2+qch2.*i;
           middlevalue2=ifft(middlevalue1, fftlen);
           ich3=real(middlevalue2);
           qch3=imag(middlevalue2);
    
    %guard interval insertion----------------------------------------------
          ich4=zeros(fftlen+gilen, nd+1);
          qch4=zeros(fftlen+gilen, nd+1);
          ich4=[ich3(fftlen-gilen+1:fftlen, :);ich3];
          qch4=[qch3(fftlen-gilen+1:fftlen, :);qch3];
    
    %attenuation calculation-----------------------------------------------
         spow=sum(ich4.^2+qch4.^2)/nd./para;
         attn=0.5*spow*sr/br*10.^(-SNR/10);
         attn=sqrt(attn);
    
    %fading channel--------------------------------------------------------
    %[ifade,qfade,ramp,rcos,rsin]=sefade(ich4,qch4,itau,dlvl1,th1,n0,itnd1,
    %nowl,length(ich4),tstp,fd,flat);
          input=reshape(ich4,1,(fftlen+gilen)*(nd+1));
          qnput=reshape(qch4,1,(fftlen+gilen)*(nd+1));
          nsamp=length(input);
          ifade=zeros(1,nsamp);
          qfade=zeros(1,nsamp);
       for k=1:now1
        % function delay--------------------------------------------------
        %function[iout,qout]=delay*(idata,qdata,nowl,idel)-----------------
        if itau~=0
            ifade(1:itau)=zeros(1,itau);
            qfade(1:itau)=zeros(1,itau);
        end
        ifade(itau+1:nsamp)=ich4(1:nsamp-itau);
        qfade(itau+1:nsamp)=qch4(1:nsamp-itau);
        %function fade-----------------------------------------------------
        %(itmp3,qtmp3,ramp,rcos,rsin)=fade(ifade,qfade,nsamp,tstp,fd,n0,---
        %itnd1,flat);------------------------------------------------------
        if fd~=0.0
            ac0=sqrt(1.0./(2.0.*(n0+1)));
            as0=sqrt(1.0./(2.0.*n0));
            ic0=itnd1;
            
            wm=2.0.*pi.*fd;
            n=4.*n0+2;
            ts=tstp;
            wmts=wm.*ts;
            paino=pi./n0;
            
            xc=zeros(1,nsamp);
            xs=zeros(1,nsamp);
            ic=[1:nsamp]+ic0;
             
           for nn=1:n0
               cwn=cos(cos(2.0.*pi.*nn./n).*ic.*wmts);
               xc=xc+cos(paino.*nn).*cwn;
               xs=xs+sin(paino.*nn).*cwn;
           end
            cwmt=sqrt(2.0).*cos(ic.*wmts);
            xc=(2.0.*xc+cwmt).*ac0;
            xs=2.0.*xs.*as0;
        
            ramp=sqrt(xc.^2+xs.^2);
            rcos=xc./ramp;
            rsin=xs./ramp;
              if flat==1
                 itmp3=sqrt(xc.^2+xs.^2).*ifade(1:nsamp);
                 qtmp3=sqrt(xc.^2+xs.^2).*qfade(1:nsamp);
              else
                 itmp3=xc.*ifade(1:nsamp)-xs.*qfade(1:nsamp);
                 qtmp3=xs.*ifade(1:nsamp)+xc.*qfade(1:nsamp);
              end
          else
            itmp3=ifade;
            qtmp3=qdata;
        end
       end
       ich5=reshape(itmp3,fftlen+gilen,nd+1);
       qch5=reshape(qtmp3,fftlen+gilen,nd+1);
       %-----------------c/i reduction-------------------------------------
       ci=10;        %carrier to interference ratio------------------------
       spowintw=sum(ich5.^2+qch5.^2)/(nd)/para;
       rint=spow/spowintw*10^(-ci/10);
       ich5=ich5.*sqrt(rint);
       qch5=qch5.*sqrt(rint);
    
    %AWGN addition---------------------------------------------------------
    ich6=zeros(fftlen+gilen, nd+1);
    qch6=zeros(fftlen+gilen, nd+1);
    rx=zeros(fftlen+gilen, nd+1);
    ich6=AWGN(ich5,SNR);
    qch6=AWGN(qch5,SNR);
    rx=ich6+qch6.*i;
    
    %guard interval removal------------------------------------------------
    ich7=zeros(fftlen, nd+1);
    qch7=zeros(fftlen, nd+1);
    ich7=[ich6(gilen+1:gilen+fftlen,:)];
    qch7=[qch6(gilen+1:gilen+fftlen,:)];
    
    %FFT-------------------------------------------------------------------
    tx1=zeros(fftlen, nd+1);
    tx2=zeros(fftlen, nd+1);
    ich8=zeros(fftlen, nd+1);
    qch8=zeros(fftlen, nd+1);
    tx1=ich7+qch7.*i;
    tx2=fft(tx1);
    ich8=real(tx2);
    qch8=imag(tx2);
    
    %--------------------fading compensation by CE symbol-----------------
    ce=1;
    ice0=ich2(:,ce);
    qce0=qch2(:,ce);
    
    ice1=ich8(:,ce);
    qce1=qch8(:,ce);
    
    iv=real((1./(ice1.^2+qce1.^2)).*(ice0+i.*qce0).*(ice1-i.*qce1));
    qv=imag((1./(ice1.^2+qce1.^2)).*(ice0+i.*qce0).*(ice1-i.*qce1));
    
    ieqv1=[iv iv iv iv iv iv iv];
    qeqv1=[qv qv qv qv qv qv qv];
    
    icompen=real((ich8+i.*qch8).*(ieqv1+i.*qeqv1));
    qcompen=imag((ich8+i.*qch8).*(ieqv1+i.*qeqv1));
    ich8=icompen;
    qch8=qcompen;
    
    %------------------CE symbol removal----------------------------------
    ich9=ich8(:,knd+1:nd+1);
    qch9=qch8(:,knd+1:nd+1);
    
    %-------------------crdemapping  DC and pilot data removal------------
    ich10(1:26,:)=ich9(2:27,:);
    qch10(1:26,:)=qch9(2:27,:);
    
    ich10(27:52,:)=ich9(39:64,:);
    qch10(27:52,:)=qch9(39:64,:);
    %QPSK demodulation-----------------------------------------------------
    outdata=zeros(para, nd*ml);
        
      outdata((1:para),(1:ml:ml*nd-1))=ich10((1:para),(1:nd))>=0;
      outdata((1:para),(2:ml:ml*nd))=qch10((1:para),(1:nd))>=0;   
      
   %parallel to serial conversion------------------------------------------
   demodata1=reshape(outdata,1,para*nd*ml);
   
   %bit error rate---------------------------------------------------------
   noe2=sum(abs(demodata1-seldata));
   nod2=length(seldata);
   noe=noe+noe2;
   nod=nod+nod2;
   if noe2~=0
       eop=eop+1;
   else
       eop=eop;
   end
   eop;
   nop=nop+1;
end
per=eop/nop;
ber=noe/nod;

BER(count,p)=ber;
p=p+1;
    end
end

average=zeros(SNR_end-SNR_start+1,1);
q=1;
for x=SNR_start:SNR_end
    average(q)=mean(BER(:,q));
    q=q+1;
end

SNR=SNR_start:SNR_end;
figure(1);
plot(SNR, average);
xlabel('SNR');
ylabel('BER -AWGN');
title('BER');
grid on;