www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/sui3/pilotLS_Lmmse_Ml.m
% 清除内存，计时开始 SUI3 信道
clear
tic;
%设置参数
totalwords=19200;
numusers=4;
wordsize=2;
linktype=0; %下行链路
ifftsize=256;
procgain=16; %扩频码的长度
guardtime=16;
guardtype=2;
delay=2;
frameguard=ifftsize+guardtime; % Guard Time between successive frames (one symbol period)
% 产生要发送的双极性二进制随机数
seed=1234;
rand('seed',seed); % Set to new seed
seqnumlist = randperm(procgain);
Datatx = zeros(numusers,totalwords);
basesignal0=zeros(1,totalwords*procgain);
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
%扩频
%ofdm调制
NumCarr =192; %取得载波数
%=============
% 每载波要传输多少数据
%=============
numsymb = ceil(length(basesignal0)/NumCarr);%numsymb=1600
%如果传输的数据不是载波的整数倍，则在后面补零
if length(basesignal0)/192 ~= numsymb,
DataPad = zeros(1,numsymb*NumCarr);
DataPad(1:length(basesignal0)) = basesignal0;
basesignal0 = DataPad;
end
clear DataPad;
%PilotIndex=[45 69 93 117 141 165 189 213]; % pilot interval=24
Data=ones(256,numsymb); % the position of data is set as 1.
Data(129,:)=0; % DC=0
Data([1:28,230:256],:)=0; % guard band=0
Data([45:24:213],:)=4; % the position of pilot is set as 4.
Data_Pilot=find(Data==4); % the pattern of pilot
Data_Pattern=find(Data==1); %the pattern of data
Data(Data_Pattern)=basesignal0; %插入数据
%插入导频数据
Wk=[1 0 1 0 0 0 1 1]'; % DownLink pilotini
pilotini=1-2*Wk; % BPSK modulation for pilot
Pilot_Data=repmat(pilotini,[1,numsymb]);
Pilot_Data=reshape(Pilot_Data,8*numsymb,1);
Data(Data_Pilot)=Pilot_Data; % pilot insertion
%==================================
%Find the time waveform using IFFT
%==================================
BaseSignal = ifft(Data); %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=[];
bertime=[];
berreal=[];
mse1=[];
mse2=[];
mse3=[];
msetime=[];
%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
L1=delay*2;
ratio=L1/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:4: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);
%提取导频数据
Rec_Pilot=Yf(Data_Pilot);
%估计导频处的信道频域响应
H_Pilot=Rec_Pilot./Pilot_Data;
H_Pilot=reshape(H_Pilot,8,numsymb);
X=17:24:185;
XI=1:201;
H_ls = INTERP1(X,H_Pilot,XI,'linear','extrap');
%时域估计
Gp=ifft(H_Pilot);
%Gp([guardtime+1:end],:)=0;
Htime=fft(Gp,201);
H_fft=fft(Gp,201);
%LMMSE
C=beita/10^(SNR/10); % constant coefficient
M=[29:229];
N=[45:24:213];
Rpp=Rhh(N,N);
Rhp=Rhh(M,N);
Wt=Rhp*inv(Rpp+C*Ip); % weight matrix of LMMSE estimation
H_lmmse=Wt*H_Pilot;
%ML
F=dftmtx(256);
L=delay*2+1;
Fh=F(N,[1:L]);
Fu=F(M,[1:L]);
H_ml=Fu*inv(Fh'*Fh)*Fh'*H_Pilot;
%计算mse
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
msetime=[msetime,mse(abs(H_real1-H_fft))]
Rx1=Yf;
Rx2=Yf;
Rx3=Yf;
Rxtime=Yf;
Rxreal=Yf;
Rx1([29:229],:)=Rx1([29:229],:)./H_ls;
Rx2([29:229],:)=Rx2([29:229],:)./H_lmmse;
Rx3([29:229],:)=Rx3([29:229],:)./H_ml;
Rxtime([29:229],:)=Rxtime([29:229],:)./Htime;
Rxreal([29:229],:)=Rxreal([29:229],:)./H_real1;
DemSig1=Rx1(Data_Pattern);
DemSig2=Rx2(Data_Pattern);
DemSig3=Rx3(Data_Pattern);
DemSigtime=Rxtime(Data_Pattern);
DemSigreal=Rxreal(Data_Pattern);
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));
DemSignaltime=reshape(DemSigtime,1,size(DemSigtime,1)*size(DemSigtime,2));
DemSignalreal=reshape(DemSigreal,1,size(DemSigreal,1)*size(DemSigreal,2));
%解扩
subsignal1=[];
subsignal2=[];
subsignal3=[];
subsignaltime=[];
subsignalreal=[];
for k=1:numusers
seqnum = seqnumlist(k);
[datarx1, subsignal10] = reccdma(DemSignal1,procgain,seqnum,linktype);
subsignal1=[subsignal1;subsignal10];
[datarx2, subsignal20] = reccdma(DemSignal2,procgain,seqnum,linktype);
subsignal2=[subsignal2;subsignal20];
[datarx3, subsignal30] = reccdma(DemSignal3,procgain,seqnum,linktype);
subsignal3=[subsignal3;subsignal30];
[datarxtime, subsignaltime0]=reccdma(DemSignaltime,procgain,seqnum,linktype);
subsignaltime=[subsignaltime;subsignaltime0];
[datarxreal, 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);
Datarxtime=invmapping(subsignaltime,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);
bertime0=err(Datatx,Datarxtime,totalwords,numusers);
berreal0=err(Datatx,Datarxreal,totalwords,numusers);
ber1=[ber1,ber10]
ber2=[ber2,ber20]
ber3=[ber3,ber30]
bertime=[bertime,bertime0]
berreal=[berreal,berreal0]
end
figure
subplot(1,2,1)
k=0:4:20
semilogy(k,ber1,'-r*')
hold on
semilogy(k,ber2,'-ro')
semilogy(k,ber3,'-rs')
semilogy(k,bertime,'-gx')
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:4:20])
hold off
subplot(1,2,2)
k=0:4:20
semilogy(k,mse1,'-r*')
hold on
semilogy(k,mse2,'-ro')
semilogy(k,mse3,'-rs')
semilogy(k,msetime,'-bx')
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:4:20])
hold off
toc