cdma.m

主要介绍了基于MATLAB环境下的CDMA无线蜂窝系统的仿真和关键物理层技术的仿真。

clear all; 
clc; 
 
%----------------------preparation part ----------------------------- 
sr=256000.0;                            % symboll rate 
ml=2;                                       % number of modulation levels 
br=sr*ml;  % bit rate 
nd=100;    % number of symbol 
ebno=8;    % Ebno 
 
%-----------------------------------------------Filter initialization----------------- 
irfn=21; 
IPOINT=8; 
alfs=0.5; 
[xh]=hrollfcoef(irfn,IPOINT,sr,alfs,1);       % Transmitting Filter 
[xh2]=hrollfcoef(irfn,IPOINT,sr,alfs,0);      % Received Filter 
 
%------------------------------ Spreading code initialization-------------------- 
user=1;   % number of users 
seq=1;  % 1:m-sequence     2:Gold    3:orthogonal Gold 
stage=3; 
ptap1=[1 3]; 
ptap2=[2 3]; 
regi1=[1 1 1]; 
regi2=[1 1 1]; 
 
% ----------------------------- Generation of the spreading code -------- 
switch seq 
    case 1 
        code=mseq(stage,ptap1,regi1,user); 
    case 2 
        m1=mseq(stage,ptap1,regi1); 
        m2=mseq(stage,ptap2,regi2); 
        code=goldseq(m1,m2,user); 
    case 3 
        m1=mseq(stage,ptap1,regi1); 
        m2=mseq(stage,ptap2,regi2); 
        code=[goldseq(m1,m2,user),zeros(user,1)]; 
end 
code=code*2-1; 
clen=length(code); 
 
% ---------------------------------- Fading initialization --------------- 
rfade=0;   % Rayleigh fading  0:nothing  1:consider 
itau=[0,8];   % delay time 
dlvl1=[0.0,40.0];  % attenuation level 
n0=[6,7]; 
th1=[0.0,0.0]; 
itnd1=[3001,4004]; 
now1=2; 
tstp=1/sr/IPOINT/clen;  % time resolution 
fd=160; 
flat=1; 
itnde1=nd*IPOINT*clen*30; 
 
% ----------------------------- Start calculation ------------- 
nloop=1000; 
noe=0; 
nod=0; 
 
for ii=1:nloop 
     
    % -------------------------- Transmitter --------------------- 
    data=rand(user,nd*ml)>0.5; 
    [ich,qch]=qpskmod(data,user,nd,ml); 
     
    [ich1,qch1]=spread(ich,qch,code); 
     
    [ich2,qch2]=compoversamp2(ich1,qch1,IPOINT); 
    [ich3,qch3]=compconv2(ich2,qch2,xh); 
     
    if user==1 
        ich4=ich3; 
        qch4=qch3; 
    else 
        ich4=sum(ich3); 
        qch4=sum(qch3); 
    end 
 
 % -------------------------- Fading channel ------------------- 
 if rfade==0 
     ich5=ich4; 
     qch5=qch4; 
 else 
     [ich5,qch5]=sefade(ich4,qch4,itau,dlvl1,th1,n0,itnd1,... 
         now1,length(ich4),tstp,fd,flat); 
     itnd1=itnd1+itnde1; 
 end 
  
 % ---------------------------------- Receiver  ----------------------------- 
 spow=sum(rot90(ich3.^2+qch3.^2))/nd; 
 attn=sqrt(0.5*spow*sr/br*10^(-ebno/10)); 
  
 [ich6,qch6]=comb2(ich5,qch5,attn); 
  
 [ich7,qch7]=compconv2(ich6,qch6,xh2); % filter 
  
 samp1=irfn*IPOINT+1; 
 ich8=ich7(:,samp1:IPOINT:IPOINT*nd*clen+samp1-1); 
 qch8=qch7(:,samp1:IPOINT:IPOINT*nd*clen+samp1-1); 
  
 [ich9 qch9]=despread(ich8,qch8,code); 
  
 demodata=qpskdemod(ich9,qch9,user,nd,ml);   % QPSK demodulation 
  
 %--------------------- BER ----------------------------- 
 noe2=sum(sum(abs(data-demodata))); 
 nod2=user*nd*ml; 
 noe=noe+noe2; 
 nod=nod+nod2; 
  
%  fprintf('%d\t%e\n',ii,noe2/nod2); 
end 
 
% ------------------------------Data file -------------------------- 
 
ber=noe/nod; 
fprintf('%d\t%d\t%d\t%e\n',ebno,noe,nod,ber); 
 
fid=fopen('BER.mat','a'); 
fprintf(fid,'%d\t%d\t%d\t%e\n',ebno,noe,nod,ber); 
 
fclose(fid); 
各子程序及其功能如下： 
1、添加高斯白噪声子程序： 
function [iout,qout]=comb2(idata,qdata,attn) 
% 
% attn: attenuation level causedby ebno or C/N 
% 
v=length(idata); 
h=length(attn); 
 
iout=zeros(h,v); 
qout=zeros(h,v); 
 
for ii=1:h 
    iout(ii,:)=idata+randn(1,v)*attn(ii); 
    qout(ii,:)=qdata+randn(1,v)*attn(ii); 
end 
2、卷积变换子程序 
function [iout,qout]=compconv2(idata,qdata,filter) 
 
% : filter tap coefficients. 
iout=conv2(idata,filter); 
qout=conv2(qdata,filter); 
 
% end of file. 
3、抽样子程序： 
function [iout,qout]=compoversamp2(iin,qin,sample) 
% 
% iin : input Ich sequence 
% qin :input Qch sequence 
% iout: ich output data sequence 
% qout : qch output data sequence 
% sample : Number of oversamples. 
 
%-------------------------------------------------------------- 
[h,v]=size(iin); 
 
iout=zeros(h,v*sample); 
qout=zeros(h,v*sample); 
 
iout(:,1:sample:1+sample*(v-1))=iin; 
qout(:,1:sample:1+sample*(v-1))=qin; 
 
% end of file. 
4、序列的自相关函数子程序： 
function [out]=autocorr(indata,tn) 
% indata: input sequence 
% tn     :  number of period 
% out   :  atuocorrelation data 
%-------------------------------------------------------------------------- 
if nargin<2 
    tn=1; 
end 
N=length(indata); 
out=zeros(1,N*tn); 
 
for ii=0:N*tn-1 
    out(ii+1)=sum(indata.*shift(indata,ii,0)); 
end 
 
% end of file 
5、序列的互相关子程序： 
function [out]=crosscorr(indata1,indata2,tn) 
% 
% indata1: input sequence1 
% indata2: input sequence2 
% tn        : number of period  
% out      : crosscorrelation data 
%---------------------------------------------------------------------- 
if nargin<3 
    tn=1; 
end  
N=length(indata1); 
out=zeros(1,N*tn); 
 
for ii=0:N*tn-1 
    out(ii+1)=sum(indata1.*shift(indata2,ii,0)); 
end 
 
% end of file. 
6、延迟子程序： 
function [out]=crosscorr(indata1,indata2,tn) 
% 
% indata1: input sequence1 
% indata2: input sequence2 
% tn        : number of period  
% out      : crosscorrelation data 
%---------------------------------------------------------------------- 
if nargin<3 
    tn=1; 
end  
N=length(indata1); 
out=zeros(1,N*tn); 
 
for ii=0:N*tn-1 
    out(ii+1)=sum(indata1.*shift(indata2,ii,0)); 
end 
 
% end of file. 
7、数据扩展子程序： 
function [iout,qout]=despread(idata,qdata,code1) 
% 
% idata : ich data sequence. 
% qdata : qch data sequence. 
% iout : ich output data sequence 
% qout : qch output data sequence 
% code1 : spread code sequence 
 
% 
switch nargin 
    case{0,1} 
        error('lack of input argument'); 
    case 2 
        code1=qdata; 
        qdata=idata; 
end 
[hn,vn]=size(idata); 
[hc,vc]=size(code1); 
 
vn=fix(vn/vc); 
iout=zeros(hc,vn); 
qout=zeros(hc,vn); 
 
for ii=1:hc 
    iout(ii,:)=rot90(flipud(rot90(reshape(idata... 
        (ii,:),vc,vn)))*rot90(code1(ii,:),3)); 
    qout(ii,:)=rot90(flipud(rot90(reshape(qdata... 
        (ii,:),vc,vn)))*rot90(code1(ii,:),3)); 
end 
 
% end of file. 
8、衰减子程序： 
function [iout,qout,ramp,rcos,rsin]=fade(idata,qdata,... 
    nsamp,tstp,fd,no,counter,flat) 
 
%****************************** varibles ********************* 
% idata: input Ich data 
% qdata: input Qch data 
% iout:  output Ich data 
% qout:  output Qch data 
% ramp:  Amplitude contaminated by fading 
% rcos:  Cosine value contaminated by fading 
% rsin:  Sine value contaminated by fading 
% nsamp: Number of samples to be simulated 
% tstp:  Minimum time resolution 
% fd:    maximum doppler frequency 
% no:    number of waves in order to generate fading 
% counter : fading counter 
% flat:     flat fading or not 
% (1-flat (only amplitude is fluctuated),0-normal(phase and amplitude are 
% fluctutated)) 
%**************************************************************** 
if fd~=0.0 
    ac0=sqrt(1.0./(2.0.*(no+1))); 
    as0=sqrt(1.0./(2.0.*no)); 
    ic0=counter; 
    pai=3.14159265; 
    wm=2.0.*pai.*fd; 
    n=4.0*no+2; 
    ts=tstp; 
    wmts=wm.*ts; 
    paino=pai./no; 
     
    xc=zeros(1,nsamp); 
    xs=zeros(1,nsamp); 
    ic=[1:nsamp]+ic0; 
     
    for nn=1:no 
        cwn=cos(cos(2.0.*pai.*nn./no).*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 
        iout=sqrt(xc.^2+xs.^2).*idata(1:nsamp); 
        qout=sqrt(xc.^2+xs.^2).*qdata(1:nsamp); 
    else 
        iout=xc.*idata(1:nsamp)-xs.*qdata(1:nsamp); 
        qout=xs.*idata(1:nsamp)+xc.*qdata(1:nsamp); 
    end 
else 
    iout=idata; 
    qout=qdata; 
end 
9、序列调制子程序： 
function [gout]=goldseq(m1,m2,n) 
% 
% m1 : m-sequence 1 
% m2:  m-sequence 2 
% n   : Number of output sequence (It can be omitted) 
% gout :output Gold sequence 
% 
if nargin<3 
    n=1; 
end  
gout=zeros(n,length(m1)); 
for ii=1:n 
    gout(ii,:)=xor(m1,m2); 
    m2=shift(m2,1,0); 
end 
 
% end of file. 
10、滤波器子程序： 
function [xh]=hrollfcoef(irfn,ipoint,sr,alfs,ncc) 
 
xi=zeros(1,irfn*ipoint+1); 
xq=zeros(1,irfn*ipoint+1); 
point=ipoint; 
tr=sr; 
tstp=1.0./tr./ipoint; 
n=ipoint.*irfn; 
mid=(n./2)+1; 
sub1=4.0.*alfs.*tr; 
 
for i=1:n 
    icon=i-mid; 
    ym=icon; 
    if icon==0.0 
        xt=(1.0-alfs+4.0.*alfs./pi).*tr; 
    else 
        sub2=16.0.*alfs.*alfs.*ym.*ym./ipoint./ipoint; 
        if sub2~=1.0 
            x1=sin(pi*(1.0-alfs)/ipoint*ym)./pi./... 
                (1.0-sub2)./ym./tstp; 
            x2=cos(pi*(1.0+alfs)/ipoint*ym)./pi.*sub1./(1.0-sub2); 
            xt=x1+x2; 
        else 
            xt=alfs.*tr.*((1.0-2.0/pi).*cos(pi/4.0/alfs)+(1.0+2.0./pi).*... 
                sin(pi/4.0/alfs))./sqrt(2.0); 
        end 
    end 
    if ncc==0 
        xh(i)=xt./ipoint./tr; 
    elseif ncc==1 
        xh(i)=xt./tr; 
    else 
        error('ncc error'); 
    end 
end 
11、m序列的调制子程序： 
% The generation function of m-sequence 
% 
% An example 
% stg=3 
% taps=[1,3] 
% inidata=[1,1,1] 
% n=2 
% 
% Programmed by yixiongshu,2006.8.9 
function [mout]=mseq(stg,taps,inidata,n) 
% stg : Number of stages 
% taps: Position of register feedback 
% inidata : Initial sequence 
% n : Number of output sequence (it can be omitted) 
% mout : output m-sequence 
% 
if nargin<4 
    n=1; 
end 
 
mout=zeros(n,2^stg-1); 
fpos=zeros(stg,1); 
fpos(taps)=1; 
 
for ii=1:2^stg-1 
    mout(1,ii)=inidata(stg); 
    num=mod(inidata*fpos,2); 
     
    inidata(2:stg)=inidata(1:stg-1); 
     
    inidata(1)=num; 
     
end 
if n>1 
    for ii=2:n 
        mout(ii,:)=shift(mout(ii-1,:),1,0); 
    end 
end 
 
% end of file. 
12、QPSK调制子程序：、 
function [demodata]=qpskdemod(idata,qdata,para,nd,ml) 
 
demodata=zeros(para,ml*nd); 
demodata((1:para),(1:ml:ml*nd-1))=idata((1:para),(1:nd))>=0; 
demodata((1:para),(2:ml:ml*nd))=qdata((1:para),(1:nd))>=0; 
function [iout,qout]=qpskmod(paradata,para,nd,ml) 
m2=ml./2; 
paradata2=paradata.*2-1; 
count2=0; 
 
for jj=1:nd 
    isi=zeros(para,1); 
    isq=zeros(para,1); 
    for ii=1:m2 
        isi=isi+2.^(m2-ii).*paradata2((1:para),ii+count2); 
        isq=isq+2.^(m2-ii).*paradata2((1:para),m2+ii+count2); 
    end 
    iout((1:para),jj)=isi; 
    qout((1:para),jj)=isq; 
    count2=count2+ml; 
end 
13、衰减子程序： 
function [iout,qout,ramp,rcos,rsin]=sefade(idata,qdata,itau,... 
    dlvl,th,n0,itn,n1,nsamp,tstp,fd,flat) 
 
%************************ variables ***************************** 
% idata input Ich data 
% qdata input  Qch data 
% iout output Ich data 
% qout output Qch data 
% ramp : Amplitude contaminated by fading 
% rcos : Cosine value contaminated by fading 
% rsin : Sine value contaminated by fading 
% itau : Delay time for each multipath fading 
% dlvl : Attenuation level for each multipath fading 
% th   : Initialized phase for each multipath fading 
% n0   : Number of waves in order to generate each multipath fading 
% itn  : Fading counter for each multipath fading 
% n1   : Number of summation for direct and delayed waves 
% nsamp: Total number of symbols  
% tstp : Minimum time resolution 
% fd   : Maximum doppler frequency 
% flat : flat fading or not 
 
iout=zeros(1,nsamp); 
qout=zeros(1,nsamp); 
 
total_attn=sum(10.^(-1.0.*dlvl./10.0)); 
for k=1:n1 
    atts=10.^(-0.05.*dlvl(k)); 
    if dlvl(k)==40.0 
        atts=0.0; 
    end 
    theta=th(k).*pi./180.0; 
    [itmp,qtmp]=delay(idata,qdata,nsamp,itau(k)); 
    [itmp3,qtmp3,ramp,rcos,rsin]=fade(itmp,qtmp,... 
        nsamp,tstp,fd,n0(k),itn(k),flat); 
    iout=iout+atts.*itmp3./sqrt(total_attn); 
    qout=qout+atts.*qtmp3./sqrt(total_attn); 
end 
14、寄存器转换子程序： 
% inregi : Vector or matrix 
% shiftr :  The amount of shift to the right. 
% shiftu : The amount of shift to the top. 
% outregi: Register output  
 
% ---------------------------------------------------- 
[h,v]=size(inregi); 
outregi=inregi; 
shiftr=rem(shiftr,v); 
shiftu=rem(shiftu,h); 
if shiftr>0 
    outregi(:,1:shiftr)=inregi(:,v-shiftr+1:v); 
    outregi(:,1+shiftr:v)=inregi(:,1:v-shiftr); 
elseif shiftr<0 
    outregi(:,1:v+shiftr)=inregi(:,1-shiftr:v); 
    outregi(:,v+shiftr+1:v)=inregi(:,1:-shiftr); 
end 
 
inregi=outregi; 
 
if shiftu>0 
    outregi(1:h-shiftu,:)=inregi(1+shiftu:h,:); 
    outregi(h-shiftu+1:h,:)=inregi(1:shiftu,:); 
elseif shiftu<0 
    outregi(1:-shiftu,:)=inregi(h+shiftu+1:h,:); 
    outregi(1-shiftu:h,:)=inregi(1:h+shiftu,:); 
end 
 
% end of file. 
15、数据扩散子程序： 
function [iout,qout]=spread(idata,qdata,code1); 
% 
% code1: spread code sequence 
% 
switch nargin 
    case{0,1} 
        error('lack of input argument'); 
    case 2 
        code1=qdata; 
        qdata=idata; 
end 
[hn,vn]=size(idata); 
[hc,vc]=size(code1); 
 
if hn>hc 
    error('lack of spread code sequences'); 
end 
 
iout=zeros(hn,vn*vc); 
qout=zeros(hn,vn*vc); 
 
for ii=1:hn                       % spread the information correctly. 
    iout(ii,:)=reshape(rot90(code1(ii,:),3)... 
        *idata(ii,:),1,vn*vc); 
    qout(ii,:)=reshape(rot90(code1(ii,:),3)... 
        *qdata(ii,:),1,vn*vc); 
end 
% end of file.