www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/sui3/preambleULLS_MMSE_ML.m
% 清除内存，计时开始 SUI3 信道，先映射后扩频
clear
tic;
%设置参数
totalwords=1920;
numusers=4;
wordsize=2;
linktype=0; %上行链路
NumCarr=256;
procgain=16; %扩频码的长度
guardtime=32;
guardtype=2;
frameguard=NumCarr+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);
H_real=zeros(NumCarr,numusers);
Rxx=[];
WW1=randsrc(256,1,[1 -1;0.5 0.5])*sqrt(2);
WW2=randsrc(256,1,[1 -1;0.5 0.5])*2;
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的扩频码
%扩频
basesignal0=tranCDMA(Datatx1,procgain,seqnum,linktype);
%=============
% 每载波要传输多少数据
%=============
Numsymb = ceil(length(basesignal0)/192);
%如果传输的数据不是数据载波的整数倍，则在后面补零
if length(basesignal0)/192 ~= Numsymb
DataPad = zeros(1,Numsymb*192);
DataPad(1:length(basesignal0)) = basesignal0;
basesignal0 = DataPad;
end
clear DataPad;
%生成数据和导频图案
Pattern=ones(NumCarr,Numsymb+2); % the position of data is set as 1.
Pattern(129,:)=0; % DC=0
Pattern([1:28,230:256],:)=0; % guard band=0
Pattern([45:24:213],:)=4;% the position of pilot is set as 4.
Pattern(:,[1,2])=0;
Data_Pattern=find(Pattern==1); %the pattern of data
%生成数据矩阵
w=round(rand(1,Numsymb)); % generate Wk
Data=zeros(size(Pattern));
for n=3:Numsymb
Data([45,93,189,213],n)=1-2*w(n); % [-84,-36,60,84]=1-2W
Data([69,117,141,165],n)=2*w(n)-1; % [-60,-12,12,36]=1-2W'
end
Data(Data_Pattern)=basesignal0; %插入用户数据
% 插入帧头
for m=1:NumCarr
p1(m,k)=WW1(m)*exp(-j*2*pi*(k-1)*guardtime*(m-1)/NumCarr);
p2(m,k)=WW2(m)*exp(-j*2*pi*(k-1)*guardtime*(m-1)/NumCarr);
end;
Data(:,1)=p1(:,k);
Data(:,2)=p2(:,k);
%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， channel have delay
%===============
% generate channel parameter which is unknown to receiver
% 抽样频率f＝4MHz，周期0.25us,信道延迟Delay=[0 0.4 0.9]us，故延迟的抽样为
% 0.4/0.25=1.6,0.9/0.25=3.6，即2和4。
Delay1=[0 0.4+0.1*(k-1) 0.9+0.4*(k-1)];
Delay(k,:)=round(Delay1/0.25);
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(k,3))];
path12=[zeros(1,Delay(k,2)) path2 zeros(1,Delay(k,3)-Delay(k,2))];
path13=[zeros(1,Delay(k,3)) path3];
RxSignal0=path11+path12+path13;
Rxx(k,:)=RxSignal0(1:length(path1));
H_real(:,k)=fft([fade(1) zeros(1,Delay(k,2)-1) fade(2) zeros(1,Delay(k,3)-Delay(k,2)-1) fade(3)].',256);
end
BER1=zeros(1,7);
BER2=zeros(1,7);
BER3=zeros(1,7);
BERreal=zeros(1,7);
MSE1=zeros(1,7);
MSE2=zeros(1,7);
MSE3=zeros(1,7);
for kk=1:numusers
ber1=[];
ber2=[];
ber3=[];
berreal=[];
mse1=[];
mse2=[];
mse3=[];
for SNR=0:5:30
RxSignal=zeros(1,size(Rxx,2));
for k=1:numusers
RxSignal=RxSignal+awgn(Rxx(k,:),10^(SNR/10),'measured',1234,'linear');
end;
%==================
% RECEIVER SECTION
%==================
% remove cyclic prefix
symbwaves=reshape(RxSignal,frameguard,Numsymb+2);
symbwaves = symbwaves(guardtime+1:frameguard,:); % Strip off the guard interval
%fft变换
Yf=fft(symbwaves);
%=======================================Estimation===============================================
%估计导频处的信道频域响应
H01=(Yf(:,2)./WW2+Yf(:,1)./WW1)/2;
%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
L=Delay(kk,3);
ratio=L/NumCarr;
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
%LS
h01=ifft(H01);
hh1=zeros(1,256);
hh1(1:guardtime)=h01((kk-1)*guardtime+1:kk*guardtime);
HH1=fft(hh1).';
%LMMSE
C=beita/10^(SNR/10); % constant coefficient
M=[1:256];
N=[1:256];
Rpp=Rhh(N,N);
Rhp=Rhh(M,N);
Wt=Rhp*inv(Rpp+C*eye(length(N))); % weight matrix of LMMSE estimation
HH2=Wt*HH1;
%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)
Nh=Delay(kk,3)+1;
Fhh=F(M,1:Nh);
Fuh=F(N,1:Nh);
HH3=Fhh*pinv(Fuh)*HH1;
H1=repmat(HH1,[1 Numsymb+2]);
H2=repmat(HH2,[1 Numsymb+2]);
H3=repmat(HH3,[1 Numsymb+2]);
H_real11=repmat(H_real(:,kk),[1 Numsymb+2]);
Hmse=H_real(:,kk);
mse1=[mse1,mse(abs(Hmse-HH1))];
mse2=[mse2,mse(abs(Hmse-HH2))];
mse3=[mse3,mse(abs(Hmse-HH3))];
Rx1=Yf./H1.*abs(H1);
Rx2=Yf./H2.*abs(H2);
Rx3=Yf./H3.*abs(H3);
Rxreal=Yf./H_real11.*abs(H_real11);
DemSig1=Rx1(Data_Pattern);
DemSig2=Rx2(Data_Pattern);
DemSig3=Rx3(Data_Pattern);
DemSigreal=Rxreal(Data_Pattern);
NN=size(DemSig1,1)*size(DemSig1,2);
DemSignal1=reshape(DemSig1,1,NN); %lS估计值
DemSignal2=reshape(DemSig2,1,NN); %MMSE估计值
DemSignal3=reshape(DemSig3,1,NN); %ML估计值
DemSignalreal=reshape(DemSigreal,1,NN);%理想估计值
%解扩
seqnum = seqnumlist(kk);
[data10,subsignal1]=recCDMA(DemSignal1,procgain,seqnum,linktype);
[data20,subsignal2]=recCDMA(DemSignal2,procgain,seqnum,linktype);
[data30,subsignal3]=recCDMA(DemSignal3,procgain,seqnum,linktype);
[datareal0,subsignalreal]=recCDMA(DemSignalreal,procgain,seqnum,linktype);
%去掉映射
Datarx1=invmapping(subsignal1,mapping,wordsize);
Datarx2=invmapping(subsignal2,mapping,wordsize);
Datarx3=invmapping(subsignal3,mapping,wordsize);
Datarxreal=invmapping(subsignalreal,mapping,wordsize);
%计算误码率
ber10=erruplink(Datatx(kk,:),Datarx1,totalwords);
ber20=erruplink(Datatx(kk,:),Datarx2,totalwords);
ber30=erruplink(Datatx(kk,:),Datarx3,totalwords);
berreal0=erruplink(Datatx(kk,:),Datarxreal,totalwords);
ber1=[ber1,ber10];
ber2=[ber2,ber20];
ber3=[ber3,ber30];
berreal=[berreal,berreal0];
end;
BER1=BER1+ber1
BER2=BER2+ber2
BER3=BER3+ber3
BERreal=BERreal+berreal
MSE1=MSE1+mse1
MSE2=MSE2+mse2
MSE3=MSE3+mse3
end;
BER1=BER1./(numusers*totalwords)
BER2=BER2./(numusers*totalwords)
BER3=BER3./(numusers*totalwords)
BERreal=BERreal./(numusers*totalwords)
MSE1=MSE1/numusers
MSE2=MSE2/numusers
MSE3=MSE3/numusers
figure
subplot(2,2,1)
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('LS','MMSE','ML','理想')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:5:30])
subplot(2,2,2)
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('LS','MMSE','ML')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:5:30])
hold off
toc