ofdm.m

OFDM在MATLAB中的仿真程序

%ofdm simulation
%-------------------------preparation part---------------------------------
para=52;
fftlen=64;
noc=53;
nd=6;

knd=1;

ml=2;
sr=250000;
br=sr*ml;
gilen=16;
ebno=1000;

%------------------------fading initialization-----------------------------
tstp=1/sr/(fftlen+gilen);
itau=[0];
dlvl1=[0];
n0=[6];
th1=[0.0];
itnd1=[1000];
now1=1;
fd=150;
flat=0;
itnd0=nd*(fftlen+gilen)*10;

%----------------------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);
%kmod=1/sqrt(2);
for iii=1:nloop
    %for i=1:para*nd*ml
     %   seldata(1,i)=round(rand);
    %end
   
    seldata=rand(1,para*nd*ml)>0.5;
    
    paradata=reshape(seldata, para, nd*ml);
    %QPSK 调制
    m2=ml./2;
    paradata2=paradata.*2-1;
    count2=0;
    for i=1:nd
        isi=zeros(para,1);
        isq=zeros(para,1);
        for j=1:m2
            isi=isi+2.^(m2-j).*paradata2((1:para),j+count2);
            isq=isq+2.^(m2-j).*paradata2((1:para),j+count2+m2);
        end
        ich0((1:para),i)=isi;
        qch0((1:para),i)=isq;
        count2=count2+ml;
    end
    %CE data 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.^(-ebno/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);
    total_attn=sum(10.^(-1.0.*dlvl1./10.0));
    for k=1:now1
        atts=10.^(-0.05.*dlvl1(k));
        if dlvl1(k)==40.0
            atts=0.0;
        end
        % 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
        itnd1=itnd1+itnd0;
        ich5=reshape(itmp3,fftlen+gilen,nd+1);
        qch5=reshape(qtmp3,fftlen+gilen,nd+1);
    end
    
    %AWGN addition
    ich6=zeros(fftlen+gilen, nd+1);
    qch6=zeros(fftlen+gilen, nd+1);
    rx=zeros(fftlen+gilen, nd+1);
    SNR=15;
    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;