www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/sui3/1users/pilotLS_Lmmse_Ml2.m
% 清除内存,计时开始 SUI3 信道 clear tic; %设置参数 totalwords=1920*5; numusers=1; wordsize=2; linktype=0; %下行链路 seqnum=3; %1~64之间取值 ifftsize=256; procgain=4; %扩频码的长度 windowtype=0; guardtime=32; guardtype=2; delay=2; frameguard=ifftsize+guardtime; % Guard Time between successive frames (one symbol period) % 产生要发送的双极性二进制随机数 rand('seed',1234); datatx = floor(rand(1,totalwords)*2); datatx = datatx-1; index = find(datatx>=0); datatx(index) = datatx(index)+1; %扩频 basesignal=tranCDMA(datatx,procgain,seqnum,linktype); index = find(basesignal<0); basesignal(index) = basesignal(index)+1; %符号转换,将双极性信号变为单极性信号 %将二进制转换为2×wordsize进制,以备映射 B=bin2high(basesignal,wordsize); %ofdm调制 NumCarr =192; %取得载波数 %============= % 每载波要传输多少数据 %============= numsymb = ceil(length(B)/NumCarr); %如果传输的数据不是载波的整数倍,则在后面补零 if length(B)/NumCarr ~= numsymb, DataPad = zeros(1,numsymb*NumCarr); DataPad(1:length(B)) = B; B = DataPad; end clear DataPad; % Mapping to the signal constellation follow mapping=get80216map(2^wordsize); for i=1:length(B), ModSignal(i)=mapping(1+B(i)); end; Xf=reshape(ModSignal,NumCarr,numsymb); %PilotIndex=[45 69 93 117 141 165 189 213]; % pilot interval=24 PilotIndex=[45:24:213]; DCIndex=129; Wk=[1 0 1 0 0 0 1 1]'; % DownLink pilotini pilotini=1-2*Wk; % BPSK modulation for pilot DataIndex=[29:229]; % the other carriers except guard carriers Xf=[zeros(28,numsymb);Xf(1:16,:);zeros(1,numsymb);Xf(17:39,:);zeros(1,numsymb);Xf(40:62,:);zeros(1,numsymb);... Xf(63:85,:);zeros(1,numsymb);Xf(86:96,:);zeros(1,numsymb);Xf(97:107,:);zeros(1,numsymb);Xf(108:130,:);zeros(1,numsymb);... Xf(131:153,:);zeros(1,numsymb);Xf(154:176,:);zeros(1,numsymb);Xf(177:192,:);zeros(27,numsymb)]; Xf(PilotIndex,:)=pilotini(:,ones(1,numsymb)); % pilot insertion Cpp=diag(Xf(PilotIndex)); % diagonal matrix composed of pilot %================================== %Find the time waveform using IFFT %================================== BaseSignal = ifft(Xf); %ifft是对列进行ifft变换。 %================================= %Add a Guard Period %================================= BaseSignal0=[BaseSignal((end-guardtime+1):end,:); BaseSignal]; BaseSignal = reshape(BaseSignal0,1,size(BaseSignal0,1)*size(BaseSignal0,2)); %先取第一列,再取第二列,…… %=============== % CHANNEL MODEL %=============== % generate channel parameter which is unknown to receiver P = [0 -5 -10]; %db P = 10.^(P/10); % calculate linear power K=[1 0 0]; s2=P./(K+1); m2=P.*(K./(K+1)); fade1=sqrt(s2)*sqrt(1/2)*(1+j); fade2=m2; fade=fade1+fade2; path1=BaseSignal*fade(1); path2=BaseSignal*fade(2); path3=BaseSignal*fade(3); path11=[path1 zeros(1,delay*2)]; path12=[zeros(1,delay) path2 zeros(1,delay)]; path13=[zeros(1,delay*2) path3]; RxSignal0=path11+path12+path13; RxSignal0=RxSignal0(1:length(path1)); H_real=fft([fade(1) zeros(1,delay-1) fade(2) zeros(1,delay-1) fade(3)].',256); H_real= H_real([29:229]); H_real1=repmat(H_real,[1 numsymb]); ber1=[]; ber2=[]; ber3=[]; berreal=[]; mse1=[]; mse2=[]; mse3=[]; %LMMSE 参数设置 switch wordsize case 2 beita=1; % beita=E{|Ck|^2}*E{1/|Ck}^2}=constant case 4 beita=17/9; % beita=17/9 for 16-QAM , and beita=1 for QPSK case 6 beita=2.6854; % beita=2.6854 for 64-QAM otherwise error('error mapping type!') end Ip=eye(8); % generate an identity matrix L=delay*2; ratio=L/ifftsize; for ri=1:256 %相关函数 for ci=1:256 if ri==ci Rhh(ri,ci)=1; else Rhh(ri,ci)=(1-exp(-j*2*pi*ratio*(ri-ci)))/(j*2*pi*ratio*(ri-ci)); % attention !!! end end end for SNR=0:2:20 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); Hp=inv(Cpp)*Yf(PilotIndex,:); % find the channel response at the pilot frequency %LMMSE C=beita/10^(SNR/10); % constant coefficient Rpp=Rhh(PilotIndex,PilotIndex); Rhp=Rhh(DataIndex,PilotIndex); Wt=Rhp*inv(Rpp+C*Ip); % weight matrix of LMMSE estimation %ML Nc=256; n=(0:Nc-1)'; k=0:Nc-1; F=1/sqrt(Nc)*exp(-j*2*pi*n*k/Nc); % generate an Nc by Nc unitary FFT matrix, satisfy:F'=inv(F) % Fuh=F(DataIndex,1:3); Nh=5; Fhh=F(DataIndex,1:Nh); Fuh=F(PilotIndex,1:Nh); DemSig1=[]; DemSig2=[]; DemSig3=[]; DemSigreal=[]; for k=1:numsymb H_ls(:,k)=interp1(PilotIndex',Hp(:,k),[DataIndex]','linear','extrap'); % 'linear' outperforms 'cubic' H_lmmse(:,k)=Wt*Hp(:,k); % LMMSE estimation H_ml(:,k)=Fhh*pinv(Fuh)*Hp(:,k); DemSig1(:,k)=Yf(DataIndex,k)./H_ls(:,k); DemSig2(:,k)=Yf(DataIndex,k)./H_lmmse(:,k); DemSig3(:,k)=Yf(DataIndex,k)./H_ml(:,k); DemSigreal(:,k)=Yf(DataIndex,k)./H_real1(:,k); end esterr1=H_ls-H_real1; % LS estimation error esterr2=H_lmmse-H_real1; % LMMSE estimation error esterr3=H_ml-H_real1; % ML estimation error mse1=[mse1,mse(abs(esterr1))] % Mean Square Error of LS estimation mse2=[mse2,mse(abs(esterr2))] % Mean Square Error of LMMSE estimation mse3=[mse3,mse(abs(esterr3))] % Mean Square Error of LMMSE estimation DemSig1([PilotIndex,DCIndex]-28,:)=[]; DemSig2([PilotIndex,DCIndex]-28,:)=[]; DemSig3([PilotIndex,DCIndex]-28,:)=[]; DemSigreal([PilotIndex,DCIndex]-28,:)=[]; DemSignal1=reshape(DemSig1,1,size(DemSig1,1)*size(DemSig1,2)); % complex signal to be de demapped DemSignal2=reshape(DemSig2,1,size(DemSig2,1)*size(DemSig2,2)); DemSignal3=reshape(DemSig3,1,size(DemSig3,1)*size(DemSig3,2)); DemSignalreal=reshape(DemSigreal,1,size(DemSigreal,1)*size(DemSigreal,2)); %去掉映射 Datarx1=invmapping(DemSignal1,mapping,NumCarr,numsymb,wordsize); Datarx2=invmapping(DemSignal2,mapping,NumCarr,numsymb,wordsize); Datarx3=invmapping(DemSignal3,mapping,NumCarr,numsymb,wordsize); Datarxreal=invmapping(DemSignalreal,mapping,NumCarr,numsymb,wordsize); %将2^wordsize进制转换为二进制, Datarx10=high2bin(Datarx1); Datarx20=high2bin(Datarx2); Datarx30=high2bin(Datarx3); Datarxreal0=high2bin(Datarxreal); %解扩 [datarx1, subsignal1] = recCDMA(Datarx10,procgain,seqnum,linktype); [datarx2, subsignal2] = recCDMA(Datarx20,procgain,seqnum,linktype); [datarx3, subsignal3] = recCDMA(Datarx30,procgain,seqnum,linktype); [datarxreal, subsignalreal] = recCDMA(Datarxreal0,procgain,seqnum,linktype); %计算误码率 ber10=err(datatx,datarx1,totalwords); ber20=err(datatx,datarx2,totalwords); ber30=err(datatx,datarx3,totalwords); berreal0=err(datatx,datarxreal,totalwords); ber1=[ber1,ber10] ber2=[ber2,ber20] ber3=[ber3,ber30] berreal=[berreal,berreal0] end figure subplot(2,2,1) k=0:2:20 semilogy(k,ber1,'-r*') hold on semilogy(k,ber2,'-ro') semilogy(k,ber3,'-rs') semilogy(k,berreal,'-.b') grid on axis([0 20 10^(-4) 0.5]) xlabel('SNR/dB'),ylabel('BER') legend('LS','LMMSE','ML','理想') set (gcf,'color',[1 1 1]) set(gca,'xtick',[0:2:20]) hold off subplot(2,2,2) k=0:2:20 semilogy(k,mse1,'-r*') hold on semilogy(k,mse2,'-ro') semilogy(k,mse3,'-rs') grid on axis([0 20 10^(-4) 1]) xlabel('SNR/dB'),ylabel('MSE') legend('LS','LMMSE','ML') set (gcf,'color',[1 1 1]) set(gca,'xtick',[0:2:20]) hold off toc