dscdma_rakematlab.m

这些代码是关于DS-CDMA 的RAKE接收机MATLAB仿真的问题（附源程序.m）已运行过。

%                          异步DS-CDMA 系统Rake接收误比特率模拟程序
K=1; %用户数，为减小仿真时间，仅考虑单个用户扩频系统。
Lc=31;%用户扩频序列周期
code_type=2;  %用户扩频序列类型；0 － 随机序列，1 － 正交序列，2 － Gold码
code_seed=rand('state'); %随机序列生成时的Seed赋值
codes=(cdma_codes(K,0,code_type,Lc,code_seed)); %调用方法见子程序。
Q=10;%扩频序列码片过抽样因子，必须满足Nyquist抽样率。
%对于矩形码片，Q应尽可能大，考虑复杂性，选择Q=10
Mpmax=5;%用户多径信道的最大路径数。
Mps=Mpmax*ones(K,1);%假设每个用户路径数相等。当然也可以设定路径数是随机数，
%如服从[1,Mpmax]间的均匀分布。

%随机产生各用户多径时延，并且保证多条路径便于分集，路径间时延间隔大于等于一个码片周期。
while (1)
    TausOri=rand(K,Mpmax)*Lc;
TausOri=sort(TausOri,2);
TausOri=floor(TausOri*Q)/Q; %考虑到时延不一定正好是抽样间隔的整数倍，这时需要内插函数
%来生成抽样时刻的信号值，比较复杂。为了减少运算量，便于编程，微调各个时延，使得
%时延正好是抽样间隔的整数倍。 
    %如果路径间时延间隔不满足条件，重新产生原始的各用户多径时延。
Tausdelta=TausOri(:,2:Mpmax)-TausOri(:,1:Mpmax-1);
    TuasSTD=sum(sum(Tausdelta>ones(K,Mpmax-1),1),2);
    if TuasSTD==K*(Mpmax-1)
        break
    end
end
TausSam=Taus*Q;%将各时延转化为需要移位的抽样间隔数
    
%多径信道模型：路径衰落功率采用等功率和指数功率模型，路径衰落相位服从[0,2π]上的均匀分布。 
%FadingsOri=kron(0:-1:1-Mpmax,ones(K,1)); %单位dB; 指数功率模型，功率为：0,－1…dB
FadingsOri=zeros(K,Mpmax); % 单位dB; 等功率模型, 功率均为0dB
Phases=rand(K,Mpmax)*2*pi;  %随机产生衰落相位
Fadings=sqrt(10.^(FadingsOri/10)).*exp(-j*Phases);%将路径衰落功率转化为衰落因子

%假设各用户路径总功率在Rake接收机处均相等，为此，对各用户发射功率比进行相应调整。
FadingRatio=abs(diag(Fadings*Fadings')/(Fadings(1,:)*Fadings(1,:)'));
NFR=1./FadingRatio; %NFR(2:K)=P2/P1=...=PK/P1   NFR(1)=P1/P1=1, 用户发射功率比值
NFRdB=10*log10(NFR); 

%用文件保存系统模拟参数
filename='CDMARakedataEqual.m';
fid=fopen(filename,'a+');
fprintf(fid,'%s ','???SimuParas: '); %’???’在文件中显示为红色，便于辨认。
SimuPara=[Lc,Taus,FadingsOri(1,:),Phases(1,:)];
for dataID=1:length(SimuPara)
    fprintf(fid,'%f  ',SimuPara(dataID));
end
fprintf(fid,'\n');
fclose(fid);

%定义和设置信噪比，信噪比定义为用户路径功率和与复加性高斯白噪声的比值
SNRdBs=-24:-14; %在相同系统参数下（包括用户数、扩频序列、多径信道等），比较不同信噪比的
%性能才有意义。
nSNRdBs=length(SNRdBs);
%逐一模拟不同信噪比下的Rake接收误比特率性能。
for SNRID=1:nSNRdBs
    SNRdB=SNRdBs(SNRID);
%信噪比定义式：SNR=Fadings(1,:)*(Fadings(1,:))'/sigma2;
    sigma2=abs(Fadings(1,:)*(Fadings(1,:))')/(10^(SNRdB/10));%CAWGN方差
    sigma=sqrt(sigma2);
    SNRbranch=abs(Fadings(1,:)).^2/sigma2;
    SNRdBbranch=10*log10(SNRbranch); %各路径信噪比
    
    nbitsAll=0; %仿真时发送的用户比特总数
    nerrorbits=zeros(Mps(1),1);%在不同Rake分集阶数下的误比特数
    times=0; %为了保证足够的误比特数，发送的用户比特总数会很大（最高达到1百万），不可能
%一次性产生所有的数据，必须逐步产生。这里以N为单位，逐步产生的次数。
  while (1) 
    N=1000;%一次性发送的用户比特数目
    bits=2*rem(unidrnd(2,K,N),2)-1; %随机产生双极性用户发送比特
    Signal=zeros(1,(N+1)*Lc*Q);%用户发送信号经过多径衰落后的信号，不包括CAWGN
    for k=1:K
        bitsChip=kron(bits(k,:),codes(k,:));%对用户比特进行扩频，比特级扩频到码片级
        ChipsSam=kron(bitsChip,ones(1,Q));%对码片进行过抽样，码片级抽样到样值级
        for pathID=1:Mps(k)
            TauSam=round(TausSam(k,pathID));
            Signal(TauSam+1:TauSam+N*Lc*Q)=Signal(TauSam+1:TauSam+N*Lc*Q)+...
                sqrt(NFR(k))*Fadings(k,pathID)*ChipsSam; %用户多径信号进行叠加
        end
    end
    %根据噪声方差，产生CAWGN
    Noise=sigma/(sqrt(2))*[1,j]*randn(2,(N+1)*Lc*Q);
    SignalRcv=Signal+Noise; %clear Signal Noise %信号叠加噪声

    %Rake接收时的各个路径分支进行相关解扩，这里假设各路径时延已经得到正确估计，
%或者说已经同步。
    TausSamEst=TausSam;
    Rakebranch=zeros(Mps(1),N); %各个路径分支相关解扩的数据
    for pathID=1:Mps(1)
        TauSam=round(TausSamEst(1,pathID));%此处程序适用于路径不是抽样间隔的整数倍的情况
        RakeChipMF=sum(reshape(SignalRcv(TauSam+1:TauSam+N*Lc*Q),Q,N*Lc),1)/Q;%矩形码片
%匹配滤波—积分清洗滤波器
        Rakebranch(pathID,:)=codes(1,:)*reshape(RakeChipMF,Lc,N)/Lc; %各个路径分支相关解扩
    end
    
%Rake 最大比率合并
    FadingsEst=Fadings(1,:); %假设信道衰落因子已经正确估计
    RakeRcv=zeros(Mps(1),N);% 不同分集阶数下的Rake接收时的解扩数据，
    for orderID=1:Mps(1)
        %按照多径时延顺序进行相应阶数下的Rake合并
RakeRcv(orderID,:)=conj(FadingsEst(1:orderID))*Rakebranch(1:orderID,:); 
    end
    %下段前面加%的程序可以画出不同Rake分集阶数下的解扩数据的星座图，用于调试程序，
%或者根据星座图中数据块分开程度，定性比较Rake分集阶数对解扩性能的影响。定量
%计算时，需要屏蔽为注释，否则会画图过多，耗尽内存资源。
%Xmax=max(max(real(RakeRcv))); Ymax=max(max(imag(RakeRcv)));
    %XYSTD=5;
    %Xmax=ceil(Xmax/XYSTD)*XYSTD;Ymax=ceil(Ymax/XYSTD)*XYSTD;
    %figure;
    %for orderID=1:Mps(1)
%   subplot(Mps(1),1, orderID);plot(real(RakeRcv(orderID,:)),imag(RakeRcv(orderID,:)),'bo');
%   axis([-Xmax,Xmax,-Ymax,Ymax]); 
%end
    %基于已仿真的数据计算误比特率。这时误比特数可能过少，需要继续发送数据。
    RakeRcvSign=sign(real(RakeRcv));
    nbitsAll=nbitsAll+N;
    nerrorbits=nerrorbits+sum(abs(sign(RakeRcvSign-kron(bits(1,:),ones(Mps(1),1)))),2);
    BER=nerrorbits/nbitsAll;
    [SNRdB,nbitsAll/N,nerrorbits'] %用于监视程序运行进展情况
    
%用文件保存系统模拟数据和结果，便于后需分析
    fid=fopen(filename,'a+');
    outdata=[SNRdB,nbitsAll,nerrorbits',BER',SNRdBbranch];
    for dataID=1:length(outdata)
         fprintf(fid,'%f  ',outdata(dataID));
    end
    fprintf(fid,'\n');
    fclose(fid);
    %判断误比特率仿真是否有效。判断标准：误比特数高于一定数目；为了防止误比特率过低，
%需要发送的用户比特总数非常巨大，对总数也进行适当约束。这里误比特标准数为100，
%用户比特总数不超过1百万，也就是说，高于10-4的误比特率仿真是有效的。
    times=times+1;
    if nerrorbits(Mps(1))>=100 | times>=1000
        break
    end
  end    
end
display('Program exit!!!'); %显示程序运行结束。这时可以打开数据保存文件，对有关数据进行处理。
%                          异步DS-CDMA 系统Rake接收误比特率模拟程序
K=1; %用户数，为减小仿真时间，仅考虑单个用户扩频系统。
Lc=31;%用户扩频序列周期
code_type=2;  %用户扩频序列类型；0 － 随机序列，1 － 正交序列，2 － Gold码
code_seed=rand('state'); %随机序列生成时的Seed赋值
codes=(cdma_codes(K,0,code_type,Lc,code_seed))'; %调用方法见子程序。
Q=10;%扩频序列码片过抽样因子，必须满足Nyquist抽样率。
%对于矩形码片，Q应尽可能大，考虑复杂性，选择Q=10
Mpmax=5;%用户多径信道的最大路径数。
Mps=Mpmax*ones(K,1);%假设每个用户路径数相等。当然也可以设定路径数是随机数，
%如服从[1,Mpmax]间的均匀分布。

%随机产生各用户多径时延，并且保证多条路径便于分集，路径间时延间隔大于等于一个码片周期。
while (1)
    TausOri=rand(K,Mpmax)*Lc;
TausOri=sort(TausOri,2);
TausOri=floor(TausOri*Q)/Q; %考虑到时延不一定正好是抽样间隔的整数倍，这时需要内插函数
%来生成抽样时刻的信号值，比较复杂。为了减少运算量，便于编程，微调各个时延，使得
%时延正好是抽样间隔的整数倍。 
    %如果路径间时延间隔不满足条件，重新产生原始的各用户多径时延。
Tausdelta=TausOri(:,2:Mpmax)-TausOri(:,1:Mpmax-1);
    TuasSTD=sum(sum(Tausdelta>ones(K,Mpmax-1),1),2);
    if TuasSTD==K*(Mpmax-1)
        break
    end
end
TausSam=Taus*Q;%将各时延转化为需要移位的抽样间隔数
    
%多径信道模型：路径衰落功率采用等功率和指数功率模型，路径衰落相位服从[0,2π]上的均匀分布。 
%FadingsOri=kron(0:-1:1-Mpmax,ones(K,1)); %单位dB; 指数功率模型，功率为：0,－1…dB
FadingsOri=zeros(K,Mpmax); % 单位dB; 等功率模型, 功率均为0dB
Phases=rand(K,Mpmax)*2*pi;  %随机产生衰落相位
Fadings=sqrt(10.^(FadingsOri/10)).*exp(-j*Phases);%将路径衰落功率转化为衰落因子

%假设各用户路径总功率在Rake接收机处均相等，为此，对各用户发射功率比进行相应调整。
FadingRatio=abs(diag(Fadings*Fadings')/(Fadings(1,:)*Fadings(1,:)'));
NFR=1./FadingRatio; %NFR(2:K)=P2/P1=...=PK/P1   NFR(1)=P1/P1=1, 用户发射功率比值
NFRdB=10*log10(NFR); 

%用文件保存系统模拟参数
filename='CDMARakedataEqual.m';
fid=fopen(filename,'a+');
fprintf(fid,'%s ','???SimuParas: '); %’???’在文件中显示为红色，便于辨认。
SimuPara=[Lc,Taus,FadingsOri(1,:),Phases(1,:)];
for dataID=1:length(SimuPara)
    fprintf(fid,'%f  ',SimuPara(dataID));
end
fprintf(fid,'\n');
fclose(fid);

%定义和设置信噪比，信噪比定义为用户路径功率和与复加性高斯白噪声的比值
SNRdBs=-24:-14; %在相同系统参数下（包括用户数、扩频序列、多径信道等），比较不同信噪比的
%性能才有意义。
nSNRdBs=length(SNRdBs);
%逐一模拟不同信噪比下的Rake接收误比特率性能。
for SNRID=1:nSNRdBs
    SNRdB=SNRdBs(SNRID);
%信噪比定义式：SNR=Fadings(1,:)*(Fadings(1,:))'/sigma2;
    sigma2=abs(Fadings(1,:)*(Fadings(1,:))')/(10^(SNRdB/10));%CAWGN方差
    sigma=sqrt(sigma2);
    SNRbranch=abs(Fadings(1,:)).^2/sigma2;
    SNRdBbranch=10*log10(SNRbranch); %各路径信噪比
    
    nbitsAll=0; %仿真时发送的用户比特总数
    nerrorbits=zeros(Mps(1),1);%在不同Rake分集阶数下的误比特数
    times=0; %为了保证足够的误比特数，发送的用户比特总数会很大（最高达到1百万），不可能
%一次性产生所有的数据，必须逐步产生。这里以N为单位，逐步产生的次数。
  while (1) 
    N=1000;%一次性发送的用户比特数目
    bits=2*rem(unidrnd(2,K,N),2)-1; %随机产生双极性用户发送比特
    Signal=zeros(1,(N+1)*Lc*Q);%用户发送信号经过多径衰落后的信号，不包括CAWGN
    for k=1:K
        bitsChip=kron(bits(k,:),codes(k,:));%对用户比特进行扩频，比特级扩频到码片级
        ChipsSam=kron(bitsChip,ones(1,Q));%对码片进行过抽样，码片级抽样到样值级
        for pathID=1:Mps(k)
            TauSam=round(TausSam(k,pathID));
            Signal(TauSam+1:TauSam+N*Lc*Q)=Signal(TauSam+1:TauSam+N*Lc*Q)+...
                sqrt(NFR(k))*Fadings(k,pathID)*ChipsSam; %用户多径信号进行叠加
        end
    end
    %根据噪声方差，产生CAWGN
    Noise=sigma/(sqrt(2))*[1,j]*randn(2,(N+1)*Lc*Q);
    SignalRcv=Signal+Noise; %clear Signal Noise %信号叠加噪声

    %Rake接收时的各个路径分支进行相关解扩，这里假设各路径时延已经得到正确估计，
%或者说已经同步。
    TausSamEst=TausSam;
    Rakebranch=zeros(Mps(1),N); %各个路径分支相关解扩的数据
    for pathID=1:Mps(1)
        TauSam=round(TausSamEst(1,pathID));%此处程序适用于路径不是抽样间隔的整数倍的情况
        RakeChipMF=sum(reshape(SignalRcv(TauSam+1:TauSam+N*Lc*Q),Q,N*Lc),1)/Q;%矩形码片
%匹配滤波—积分清洗滤波器
        Rakebranch(pathID,:)=codes(1,:)*reshape(RakeChipMF,Lc,N)/Lc; %各个路径分支相关解扩
    end
    
%Rake 最大比率合并
    FadingsEst=Fadings(1,:); %假设信道衰落因子已经正确估计
    RakeRcv=zeros(Mps(1),N);% 不同分集阶数下的Rake接收时的解扩数据，
    for orderID=1:Mps(1)
        %按照多径时延顺序进行相应阶数下的Rake合并
RakeRcv(orderID,:)=conj(FadingsEst(1:orderID))*Rakebranch(1:orderID,:); 
    end
    %下段前面加%的程序可以画出不同Rake分集阶数下的解扩数据的星座图，用于调试程序，
%或者根据星座图中数据块分开程度，定性比较Rake分集阶数对解扩性能的影响。定量
%计算时，需要屏蔽为注释，否则会画图过多，耗尽内存资源。
%Xmax=max(max(real(RakeRcv))); Ymax=max(max(imag(RakeRcv)));
    %XYSTD=5;
    %Xmax=ceil(Xmax/XYSTD)*XYSTD;Ymax=ceil(Ymax/XYSTD)*XYSTD;
    %figure;
    %for orderID=1:Mps(1)
%   subplot(Mps(1),1, orderID);plot(real(RakeRcv(orderID,:)),imag(RakeRcv(orderID,:)),'bo');
%   axis([-Xmax,Xmax,-Ymax,Ymax]); 
%end
    %基于已仿真的数据计算误比特率。这时误比特数可能过少，需要继续发送数据。
    RakeRcvSign=sign(real(RakeRcv));
    nbitsAll=nbitsAll+N;
    nerrorbits=nerrorbits+sum(abs(sign(RakeRcvSign-kron(bits(1,:),ones(Mps(1),1)))),2);
    BER=nerrorbits/nbitsAll;
    [SNRdB,nbitsAll/N,nerrorbits'] %用于监视程序运行进展情况
    
%用文件保存系统模拟数据和结果，便于后需分析
    fid=fopen(filename,'a+');
    outdata=[SNRdB,nbitsAll,nerrorbits',BER',SNRdBbranch];
    for dataID=1:length(outdata)
         fprintf(fid,'%f  ',outdata(dataID));
    end
    fprintf(fid,'\n');
    fclose(fid);
    %判断误比特率仿真是否有效。判断标准：误比特数高于一定数目；为了防止误比特率过低，
%需要发送的用户比特总数非常巨大，对总数也进行适当约束。这里误比特标准数为100，
%用户比特总数不超过1百万，也就是说，高于10-4的误比特率仿真是有效的。
    times=times+1;
    if nerrorbits(Mps(1))>=100 | times>=1000
        break
    end
  end    
end
display('Program exit!!!'); %显示程序运行结束。这时可以打开数据保存文件，对有关数据进行处理。
%                          异步DS-CDMA 系统Rake接收误比特率模拟程序