www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/sui3/1/mccdma1.m
% 清除内存，计时开始 变码长 SUI3 信道
clear
tic;
%设置参数
totalwords=1920;
wordsize=2;
linktype=0; %下行链路
ifftsize=256;
numusers=4;
procgain=16; %扩频码的长度
guardtime=16;
guardtype=2;
frameguard=ifftsize+guardtime; % Guard Time between successive frames (one symbol period)
% 产生要发送的双极性二进制随机数
ber=[];
mse_ls=[];
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);
%扩频
basesignal_0=tranCDMA(Datatx1,procgain,seqnum,linktype);
basesignal0=basesignal0+basesignal_0; %各个用户的扩频信号进行合并
end
%ofdm调制
NumCarr =192; %取得载波数
%=============
% 每载波要传输多少数据
%=============
numsymb = ceil(length(basesignal0)/NumCarr);
%如果传输的数据不是载波的整数倍，则在后面补零
if length(basesignal0)/NumCarr ~= numsymb,
DataPad = zeros(1,numsymb*NumCarr);
DataPad(1:length(basesignal0)) = basesignal0;
basesignal0 = DataPad;
end
clear DataPad;
Xf=reshape(basesignal0,NumCarr,length(basesignal0)/NumCarr);
%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
%=================================
if guardtype~=0
if guardtype == 1 %if guard period is just a zero transmission
BaseSignal=[zeros(guardtime,numsymb);BaseSignal];
elseif guardtype == 2
EndSignal = size(BaseSignal,1); %Find the number of columns in the BaseSignal
%将后guardtime拷贝，作为循环前缀，即：L＝guardtime，N＝EndSignal。
BaseSignal=[BaseSignal((EndSignal-guardtime+1):EndSignal,:); BaseSignal];
elseif guardtype == 3
EndSignal = size(BaseSignal,1); %Find the number of columns in the BaseSignal
BaseSignal=[zeros(guardtime/2,numsymb); ...
BaseSignal((EndSignal-guardtime/2+1):EndSignal,:); BaseSignal];
end
end
%BaseSignal = reshape(BaseSignal,1,size(BaseSignal,1)*size(BaseSignal,2)) %先取第一列，再取第二列，……
%===============
% CHANNEL MODEL
%===============
% generate channel parameter which is unknown to receiver
P=[0 -5 -10]; % power of each tap in dB
K=[1 0 0]; % Ricean K-factor in linear scale
N=10000; % number of independent random realizations
Fnorm=-1.5113; % gain normalization factor in dB
dop=[0.4 0.4 0.4]; % Doppler maximal frequency parameter in Hz
M=256; % number of taps of the Doppler filter
P=10.^(P/10); % linear power:[1.0000 0.3162 0.1000]
s2=P./(K+1); % variance: [0.5000 0.3162 0.1000]
m2=P.*(K./(K+1)); % constant power: [0.5000 0 0]. m2=P-s2 (direct power=total power-diffuse power)
m=sqrt(m2); % constant part
L=length(P);
paths_r=sqrt(1/2)*(randn(L,N)+j*randn(L,N)).*((sqrt(s2))'*ones(1,N));
paths_c=m'*ones(1,N);
for p=1:L
D=dop(p)/max(dop)/2; % normalize to highest Doppler
f0=[0:M*D]/(M*D); % freqency vector
psd=0.785*f0.^4-1.72*f0.^2+1.0; % psd approximation
filt=[psd(1:end-1) zeros(1,M-2*M*D) psd(end:-1:2)]; % S(f)
filt=sqrt(filt); % form S(f) to |H(f)|
filt=ifftshift(ifft(filt)); % get impulse response
filt=real(filt); % want a real-valued filter
filt=filt/sqrt(sum(filt.^2)); % normalize filter
path=fftfilt(filt,[paths_r(p,:) zeros(1,M)]);
paths_r(p,:)=path(1+M/2:end-M/2);
end;
paths=paths_r+paths_c;
paths=paths*10^(Fnorm/20);
path1=ones(frameguard,1)*paths(1,[101:numsymb+100]).*BaseSignal;
path2=ones(frameguard,1)*paths(2,[101:numsymb+100]).*BaseSignal;
path3=ones(frameguard,1)*paths(3,[101:numsymb+100]).*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));
delay=2; % multipath delay(in sample)
pathd1=[path01 zeros(1,delay*2)];
pathd2=[zeros(1,delay) path02 zeros(1,delay)];
pathd3=[zeros(1,delay*2) path03];
y=pathd1+pathd2+pathd3;
y=y(1:length(path01)); % sync in the path1
for SNR=0:4:40
RxSignal=awgn(y,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
DemSig=[];
for k=1:numsymb
H_ls(:,k)=interp1(PilotIndex',Hp(:,k),[DataIndex]','linear','extrap');
H_ideal(:,k)=fft([paths(1,k+100) zeros(1,delay-1) paths(2,k+100) zeros(1,delay-1) paths(3,k+100)].',256);
DemSig(:,k)=Yf(DataIndex,k)./H_ls(:,k); % extract effective data
end
esterr=H_ls-H_ideal(DataIndex,:); % LS estimation error
mse_ls=[mse_ls,mse(abs(esterr))]; % Mean Square Error of LS estimation
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
%解扩
subsignal=[];
for k=1:numusers
seqnum = seqnumlist(k);
[data10,subsignal10]=recCDMA(DemSig,procgain,seqnum,linktype);
subsignal=[subsignal;subsignal10];
end
%去映射
Datarx=invmapping(subsignal,mapping,wordsize);
%计算误码率
ber10=err(Datatx,Datarx,totalwords,numusers);
ber=[ber,ber10];
end
figure
subplot(2,2,1)
k=0:4:40
semilogy(k,ber,'-b*','MarkerSize',5)
grid on
axis([0 40 10^(-4) 1])
xlabel('SNR/dB'),ylabel('BER')
legend('numusers=4')
set (gcf,'color',[1 1 1])
subplot(2,2,2)
k=0:4:40
semilogy(k,mse_ls,'-b*','MarkerSize',5)
grid on
axis([0 40 10^(-4) 1])
xlabel('SNR/dB'),ylabel('mse')
legend('numusers=4')
set (gcf,'color',[1 1 1])
toc