www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/sui3/1users/pilotLS.m
% 清除内存，计时开始 SUI3 信道
clear
tic;
%设置参数
totalwords=1920*6;
numusers=4;
wordsize=2; %QPSK
linktype=0; %上行链路
seqnum=3; %1～64之间取值
ifftsize=256;
procgain=16; %扩频码的长度
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);%0和1
datatx = datatx-1;%-1和0
index = find(datatx>=0);
datatx(index) = datatx(index)+1;%为0的变为+1
ber=[];
MSE=[];
%扩频
for i=2:1:4
procgain=2^i;
basesignal=trancdma(datatx,procgain,seqnum,linktype);%用户数据与扩频码相乘 %1920*6*2*16
index = find(basesignal<0);
basesignal(index) = basesignal(index)+1; %符号转换，将双极性信号变为单极性信号
%将二进制转换为2×wordsize进制，以备映射
B=bin2high(basesignal,wordsize);%ofdm调制 %B为1行，size(basesignal,2)/wordsize列 %1920*6*2*16列
NumCarr =192; %取得载波数
%=============
% 每载波要传输多少数据
%=============
numsymb = ceil(length(B)/NumCarr);%120*16
%如果传输的数据不是载波的整数倍，则在后面补零
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]);
for SNR=0:2:24
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
for n=0:255
k=0:255;
F(n+1,:)=1/sqrt(256)*exp(-j*2*pi*n*k/256); % N by N unitary FFT matrix, note:F'=inv(F)
end
Fuh=F(DataIndex,1:3);
DemSig=[];
for k=1:numsymb
H_ls(:,k)=interp1(PilotIndex',Hp(:,k),[DataIndex]','linear','extrap');
DemSig(:,k)=Yf(DataIndex,k)./H_ls(:,k); % extract effective data
end
esterr=H_ls-H_real1; % LS estimation error
mse_ls=mse(abs(esterr)); % Mean Square Error of LS estimation
MSE=[MSE,mse_ls]
fprintf('SNR=%idB, MSE of LS estimation=%i\n',SNR,mse_ls);
DemSig([DCIndex,PilotIndex]-28,:)=[]; % DC and pilots removal
DemSig=reshape(DemSig,1,size(DemSig,1)*size(DemSig,2)); % complex signal to be de demapped
%去映射
Datarx=zeros(1,NumCarr*numsymb);
M=2^wordsize;
for i=1:NumCarr*numsymb
for j=1:M,
metrics(j)=(real(DemSig(i))-real(mapping(j)))^2+(imag(DemSig(i))-imag(mapping(j)))^2;
end;
[min_metrics decis]=min(metrics); % [a b]=min(metrics) a:minimum value，b: responding position
Datarx(i)=decis-1;
end;
%解扩
%将2×wordsize进制二进制转换为二进制，以备映射
d=Datarx;
C=dec2bin(d);%十进制转为二进制
C=C';
Datarx=reshape(C,1,2*size(C,2));
index = find(Datarx<=0);
Datarx(index) = Datarx(index)-1;%符号转换，将单性信号变为双极性信号
[datarx, subsignal] = reccdma(Datarx,procgain,seqnum,linktype);
errs=totalwords;
for i=1:totalwords,
if datatx(i)==datarx(i),
errs=errs-1;
end;
end;
ber=[ber,errs/totalwords]
end;%for SNR=0:2:24
end;%for i=2:1:4
figure
subplot(2,2,1)
k=0:2:24
semilogy(k,ber(1,1:13),'-r*')
hold on
semilogy(k,ber(1,14:26),'-bo')
semilogy(k,ber(1,27:39),'-rs')
grid on
axis([0 24 10^(-4) 0.5])
xlabel('SNR/dB'),ylabel('BER')
legend('procgain=4','procgain=8','procgain=16')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:2:24])
subplot(2,2,2)
k=0:2:24
semilogy(k,MSE(1,1:13),'-r*')
hold on
semilogy(k,MSE(1,14:26),'-bo')
semilogy(k,MSE(1,27:39),'-rs')
grid on
axis([0 24 10^(-4) 1])
xlabel('SNR/dB'),ylabel('MSE')
legend('procgain=4','procgain=8','procgain=16')
set (gcf,'color',[1 1 1])
set(gca,'xtick',[0:2:24])
toc