www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/rayleigh(移动信道)/twoDft.m

    % 清除内存，计时开始  4用户瑞利信道256点
clear 
tic;
%设置参数
numusers=4;
wordsize=2;
linktype=0;   %下行链路
NumCarr=256;
procgain=16;  %扩频码的长度
guardtime=40;
guardtype=2;
frameguard=NumCarr+guardtime;  % Guard Time between successive frames (one symbol period)
Numsymb=1600;
fc=3.5e9;  
V=75;         %V km/h
fdmax=V*fc/3e8/3.6   %fmax=V*fc/C

%ofdm调制
%=============
% 每载波要传输多少数据
%=============

%=================
% pilot settings
%=================
fnum=40;
fspace=5;
tspace=2;
%生成数据和导频图案
Pattern=ones(NumCarr,Numsymb);  % the position of data is set as 1.
Pattern(129,:)=0;                 % DC=0
Pattern([1:28,230:256],:)=0;      % guard band=0
Pattern([31:fspace:226],[1:tspace:Numsymb])=4;       % the position of pilot is set as 4.
Pattern_Pilot=find(Pattern==4);   % the pattern of pilot
Data_Pattern=find(Pattern==1);   %the pattern of data
Data=zeros(size(Pattern));
% 产生要发送的随机数
seed=1234;
rand('seed',seed);	                    	% Set to new seed
seqnumlist = randperm(procgain);
totalwords=length(Data_Pattern)/procgain;
Datatx = zeros(numusers,totalwords);
basesignal0=zeros(1,length(Data_Pattern));
for k = 1:numusers,  
    Datatx(k,:)=floor(rand(1,totalwords)*2^wordsize);%产生2^wordsize进制发送数据
    %进行映射
     mapping=get80216map(2^wordsize);
     B=Datatx(k,:);
     for i=1:length(B),
     Datatx1(i)=mapping(1+B(i));
     end;
    seqnum=seqnumlist(k);%取得用户k的扩频码
%扩频
    basesignal_0=tranCDMA(Datatx1,procgain,seqnum,linktype);
    basesignal0=basesignal0+basesignal_0; %各个用户的扩频信号进行合并
 end

Data(Data_Pattern)=basesignal0;   %插入用户数据
%生成导频数据
Pilot_Data=randsrc(length(Pattern_Pilot),1,[1 -1;0.5 0.5])*sqrt(2);
Data(Pattern_Pilot)=Pilot_Data;
%==================================
%Find the time waveform using IFFT
%==================================
BaseSignal = ifft(Data);       %ifft是对列进行ifft变换。
%=================================
%Add a Guard Period
%=================================
BaseSignal=[BaseSignal((end-guardtime+1):end,:); BaseSignal];
%BaseSignal = reshape(BaseSignal0,1,size(BaseSignal0,1)*size(BaseSignal0,2));   %先取第一列，再取第二列，……

%===============
% CHANNEL MODEL， channel have delay
%===============
fade=Rayleigh(fdmax);
path1=ones(frameguard,1)*fade(1,[1+5000:Numsymb+5000]).*BaseSignal;
path2=ones(frameguard,1)*fade(2,[1+5000:Numsymb+5000]).*BaseSignal;
path3=ones(frameguard,1)*fade(3,[1+5000:Numsymb+5000]).*BaseSignal;
path4=ones(frameguard,1)*fade(4,[1+5000:Numsymb+5000]).*BaseSignal;
path5=ones(frameguard,1)*fade(5,[1+5000:Numsymb+5000]).*BaseSignal;
path6=ones(frameguard,1)*fade(6,[1+5000:Numsymb+5000]).*BaseSignal;
path01=reshape(path1,1,size(path1,1)*size(path1,2));
path02=reshape(path2,1,size(path2,1)*size(path2,2));
path03=reshape(path3,1,size(path3,1)*size(path3,2));
path04=reshape(path4,1,size(path4,1)*size(path4,2));
path05=reshape(path5,1,size(path5,1)*size(path5,2));
path06=reshape(path6,1,size(path6,1)*size(path6,2));
%the delay is [0 0.31 0.71 1.09 1.73 2.51] us 抽样频率为4MHz,[0 1.24 2.84 4.36 6.92 10]samples
path11=[path01 zeros(1,10)];    % the largest delay is 10 sample                                                                                                                          
path12=[zeros(1,1) path02 zeros(1,9)];
path13=[zeros(1,3) path03 zeros(1,7) ];
path14=[zeros(1,4) path04 zeros(1,6)];
path15=[zeros(1,7) path05 zeros(1,3)];
path16=[zeros(1,10) path06];
RxSignal0=path11+path12+path13+path14+path15+path16;
RxSignal0=RxSignal0(1:length(path01));

H_real=zeros(NumCarr,Numsymb);      % the real channel impulse response
for k=1:NumCarr
    H_real(k,:)=H_real(k,:)+fade(1,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*0/NumCarr)+fade(2,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*1/NumCarr)+...
                            fade(3,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*3/NumCarr)+fade(4,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*4/NumCarr)+...
                            fade(5,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*7/NumCarr)+fade(6,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*10/NumCarr);
end
H_real= H_real([29:229],:);
ber1=[];
%ber2=[];
%ber3=[];
berreal=[];
mse1=[];
%mse2=[];
%mse3=[];
for SNR=0:5:30
RxSignal=awgn(RxSignal0,10^(SNR/10),'measured',1234,'linear');
%==================
% RECEIVER SECTION
%==================
% remove cyclic prefix
symbwaves=reshape(RxSignal,frameguard,Numsymb);
symbwaves = symbwaves(guardtime+1:frameguard,:); % Strip off the guard interval
%fft变换
Yf=fft(symbwaves);
%=======================================Estimation===============================================
%提取导频数据
	Rec_Pilot=Yf(Pattern_Pilot);
%估计导频处的信道频域响应
	H_Pilot=Rec_Pilot./Pilot_Data;
    H_Pilot=reshape(H_Pilot,fnum,length(H_Pilot)/fnum);
for m=1:size(H_Pilot,2);
     H_fdir(:,m)=interp1([3:fspace:198]',H_Pilot(:,m),[1:201]','linear','extrap');   % linear interpolation in frequency domain
end
for k=1:201
    H1(k,:)=interp1([1:tspace:Numsymb],H_fdir(k,:),[1:Numsymb],'linear','extrap'); % linear interpolation in time domain
end
   
 mse1=[mse1,mse(abs(H_real-H1))]
% mse2=[mse2,mse(abs(H_real-H2))]
 %mse3=[mse3,mse(abs(H_real-H3))]
  H10=abs(H1);
   Rx1=Yf;
   %Rx2=Yf;
   %Rx3=Yf;
   Rxreal=Yf;
   Rx1([29:229],:)=Rx1([29:229],:)./H1.*H10;
   %Rx2([29:229],:)=Rx2([29:229],:)./H2;
   %Rx3([29:229],:)=Rx3([29:229],:)./H3;
   Rxreal([29:229],:)=Rxreal([29:229],:)./H_real.*abs(H_real);
    
   DemSig1=Rx1(Data_Pattern);
   %DemSig2=Rx2(Data_Pattern);
   %DemSig3=Rx3(Data_Pattern);
   DemSigreal=Rxreal(Data_Pattern);
   
DemSignal1=reshape(DemSig1,1,size(DemSig1,1)*size(DemSig1,2));   %LS估计值
%DemSignal2=reshape(DemSig2,1,size(DemSig2,1)*size(DemSig2,2));   %MMSE估计值
%DemSignal3=reshape(DemSig3,1,size(DemSig3,1)*size(DemSig3,2));   %ML估计值
DemSignalreal=reshape(DemSigreal,1,size(DemSigreal,1)*size(DemSigreal,2));%理想估计值
%解扩
subsignal1=[];
%subsignal2=[];
%subsignal3=[];
subsignalreal=[];
for k=1:numusers
    seqnum = seqnumlist(k);
[data10,subsignal10]=recCDMA(DemSignal1,procgain,seqnum,linktype);
subsignal1=[subsignal1;subsignal10];
%[data20,subsignal20]=recCDMA(DemSignal2,procgain,seqnum,linktype);
%subsignal2=[subsignal2;subsignal20];
%[data30,subsignal30]=recCDMA(DemSignal3,procgain,seqnum,linktype);
%subsignal3=[subsignal3;subsignal30];
[datareal0,subsignalreal0]=recCDMA(DemSignalreal,procgain,seqnum,linktype);
subsignalreal=[subsignalreal;subsignalreal0];
end
%去掉映射
Datarx1=invmapping(subsignal1,mapping,wordsize);
%Datarx2=invmapping(subsignal2,mapping,wordsize);
%Datarx3=invmapping(subsignal3,mapping,wordsize);
Datarxreal=invmapping(subsignalreal,mapping,wordsize);
%计算误码率
ber10=err(Datatx,Datarx1,totalwords,numusers);
%ber20=err(Datatx,Datarx2,totalwords,numusers);
%ber30=err(Datatx,Datarx3,totalwords,numusers);
berreal0=err(Datatx,Datarxreal,totalwords,numusers)
ber1=[ber1,ber10]
%ber2=[ber2,ber20]
%ber3=[ber3,ber30]
berreal=[berreal,berreal0]
end;

figure
subplot(2,2,3)
k=0:5:30
semilogy(k,ber1,'-rh')
hold on
%semilogy(k,ber2,'-ro')
%semilogy(k,ber3,'-rx')
semilogy(k,berreal,'-.b')
grid on
axis([0 30 10^(-4) 1])
xlabel('SNR/dB'),ylabel('BER')
legend('150km/h','75km/h','3km/h')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:5:30])

subplot(2,2,4)
k=0:5:30
semilogy(k,mse1,'-rh')
hold on
%semilogy(k,mse2,'-ro')
%semilogy(k,mse3,'-rx')
grid on
axis([0 30 10^(-4) 1])
xlabel('SNR/dB'),ylabel('MSE')
legend('150km/h','75km/h','3km/h')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:5:30])

toc