www.gusucode.com > MC-CDMA系统的仿真matlab源码程序 > mc-cdma/rayleigh(移动信道)/poly_wenxian.m
% 清除内存，计时开始 6径Rayleigh信道,导频数不能取太小，取20效果较好
clear
tic;
%设置参数
order=3;
MAXSIZE=8;
totalwords=90;
numusers=4;
wordsize=2;
linktype=0; %下行链路
ifftsize=256;
procgain=16; %扩频码的长度
guardtime=32;
guardtype=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 =180; %数据载波数
%=============
% 每载波要传输多少数据
%=============
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;
%PilotIndex=[45 69 93 117 141 165 189 213]; % pilot interval=24
%生成数据和导频图案
Pilot_number=20;
Pattern=ones(ifftsize,Numsymb); % the position of data is set as 1.
Pattern(129,:)=0; % DC=0
Pattern([1:28,230:256],:)=0; % guard band=0
Pattern([34:10:224],:)=4; % the position of pilot is set as 4.
PilotIndex=find(Pattern==4); % the pattern of pilot
DataIndex=find(Pattern==1); %the pattern of data
%生成导频数据
pilotini=randsrc(length(PilotIndex),1,[1 -1;0.5 0.5])*sqrt(2);
%生成数据矩阵
Data=zeros(size(Pattern));
Data(PilotIndex)=pilotini;
Data(DataIndex)=basesignal0;
%==================================
%Find the time waveform using IFFT
%==================================
BaseSignal = ifft(Data); %ifft是对列进行ifft变换。
%将后guardtime拷贝，作为循环前缀，即：L＝guardtime，N＝EndSignal。
BaseSignal=[BaseSignal((end-guardtime+1):end,:); BaseSignal];
%BaseSignal = reshape(BaseSignal,1,size(BaseSignal,1)*size(BaseSignal,2)) %先取第一列，再取第二列，……
%===============
% CHANNEL MODEL
%===============
%fc=3.5e9;
%V=75; %V km/h
fdmax=200;%V*fc/3e8/3.6 %fmax=V*fc/C
% channel have delay
fade=Rayleigh(fdmax);
% fade=Rayleigh_CH2(fdmax);
path1=ones(frameguard,1)*fade(1,[1+5000:Numsymb+5000]).*BaseSignal;
path2=ones(frameguard,1)*fade(2,[1+5000:Numsymb+5000]).*BaseSignal;
path3=ones(frameguard,1)*fade(3,[1+5000:Numsymb+5000]).*BaseSignal;
path4=ones(frameguard,1)*fade(4,[1+5000:Numsymb+5000]).*BaseSignal;
path5=ones(frameguard,1)*fade(5,[1+5000:Numsymb+5000]).*BaseSignal;
path6=ones(frameguard,1)*fade(6,[1+5000:Numsymb+5000]).*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));
path04=reshape(path4,1,size(path4,1)*size(path4,2));
path05=reshape(path5,1,size(path5,1)*size(path5,2));
path06=reshape(path6,1,size(path6,1)*size(path6,2));
%the delay is [0 0.31 0.71 1.09 1.73 2.51] us
path11=[path01 zeros(1,10)]; % the largest delay is 10 sample
path12=[zeros(1,1) path02 zeros(1,9)];
path13=[zeros(1,3) path03 zeros(1,7) ];
path14=[zeros(1,4) path04 zeros(1,6)];
path15=[zeros(1,7) path05 zeros(1,3)];
path16=[zeros(1,10) path06];
RxSignal0=path11+path12+path13+path14+path15+path16;
RxSignal0=RxSignal0(1:length(path01));
H_real=zeros(ifftsize,Numsymb); % the real channel impulse response
for k=1:ifftsize
H_real(k,:)=H_real(k,:)+fade(1,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*0/ifftsize)+fade(2,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*1/ifftsize)+...
fade(3,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*3/ifftsize)+fade(4,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*4/ifftsize)+...
fade(5,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*7/ifftsize)+fade(6,[1+5000:Numsymb+5000])*exp(-j*2*pi*(k-1)*10/ifftsize);
end;
H_real1= H_real([29:229],:);
H_real2=H_real(PilotIndex);
h_real=ifft(H_real2);
h_real=reshape(h_real,Pilot_number,Numsymb);
ber1=[];
ber2=[];
ber3=[];
berreal=[];
mse1=[];
mse2=[];
mse3=[];
for SNR=0:5:30
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);
%=======================================Estimation===============================================
%提取导频数据
Rec_Pilot=Yf(PilotIndex);
%估计导频处的信道频域响应
H_Pilot=Rec_Pilot./pilotini;
H_Pilot=reshape(H_Pilot,Pilot_number,Numsymb);%20*1600
X=6:10:196;
XI=1:201;
H_ls=INTERP1(X,H_Pilot,XI,'linear','extrap');%201*1600
h_ls=ifft(H_Pilot);
h_p=h_ls;%20*1600
h_p1=abs(h_p).^2;
Pw=(sum(h_p1'))'/Numsymb;%20*1
[A1,I1]=sort(Pw);
h_p(I1(1:Pilot_number-MAXSIZE),:)=0;%h_MST
h_MST=h_p;
H_MST=fft(h_p,201);
%h_MST=ifft(H_MST);
C=h_p(I1(Pilot_number-MAXSIZE+1:Pilot_number),:);%6*1600
B=ones(order+1,Numsymb);
B(2,:)=(0:Numsymb-1);
if order>1
for kk=2:order
B(kk+1,:)=B(2,:).^kk;
end
end
A=C*B'*inv(B*B');
for ii=0:Numsymb-1
for jj=1:size(A,1)
switch order
case 1
h_poly(I1((Pilot_number-MAXSIZE+1)+(jj-1)),ii+1)=A(jj,2)*(ii^1)+A(jj,1);
case 2
h_poly(I1((Pilot_number-MAXSIZE+1)+(jj-1)),ii+1)=A(jj,3)*(ii^2)+A(jj,2)*(ii^1)+A(jj,1);
case 3
h_poly(I1((Pilot_number-MAXSIZE+1)+(jj-1)),ii+1)=A(jj,4)*(ii^3)+A(jj,3)*(ii^2)+A(jj,2)*(ii^1)+A(jj,1);
case 4
h_poly(I1((Pilot_number-MAXSIZE+1)+(jj-1)),ii+1)=A(jj,5)*(ii^4)+A(jj,4)*(ii^3)+A(jj,3)*(ii^2)+A(jj,2)*(ii^1)+A(jj,1);
end
end
end
h_poly(I1(1:Pilot_number-MAXSIZE),:)=0;
H_poly=fft(h_poly,201);
%h_poly=ifft(H_poly);
mse1=[mse1,mse(abs(h_real-h_ls))];
mse2=[mse2,mse(abs(h_real-h_poly))];
mse3=[mse3,mse(abs(h_real-h_MST))];
Rx1=Yf;
Rx2=Yf;
Rx3=Yf;
Rxreal=Yf;
Rx1([29:229],:)=Rx1([29:229],:)./H_ls.*abs(H_ls);
Rx2([29:229],:)=Rx2([29:229],:)./H_poly.*abs(H_poly);%最大合并比，×第i个子载波处的随机幅度
Rx3([29:229],:)=Rx3([29:229],:)./H_MST.*abs(H_MST);
Rxreal([29:229],:)=Rxreal([29:229],:)./H_real1.*abs(H_real1);
% extract data
DemSig1=Rx1(DataIndex);
DemSig2=Rx2(DataIndex);
DemSig3=Rx3(DataIndex);
DemSigreal=Rxreal(DataIndex);
DemSignal1=reshape(DemSig1,1,size(DemSig1,1)*size(DemSig1,2));
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));%理想估计值
%解扩
subsignal1=[];
subsignal2=[];
subsignal3=[];
subsignalreal=[];
for k=1:numusers
seqnum = seqnumlist(k);
[data1,subsignal10]=recCDMA(DemSignal1,procgain,seqnum,linktype);
subsignal1=[subsignal1;subsignal10];
[data2,subsignal20]=recCDMA(DemSignal2,procgain,seqnum,linktype);
subsignal2=[subsignal2;subsignal20];
[data3,subsignal30]=recCDMA(DemSignal3,procgain,seqnum,linktype);
subsignal3=[subsignal3;subsignal30];
[datareal0,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);
Datarxreal=invmapping(subsignalreal,mapping,wordsize);
%计算误码率
ber10=err(Datatx,Datarx1,totalwords,numusers);
ber20=err(Datatx,Datarx2,totalwords,numusers);
ber30=err(Datatx,Datarx3,totalwords,numusers);
berreal0=err(Datatx,Datarxreal,totalwords,numusers);
ber1=[ber1,ber10];
ber2=[ber2,ber20];
ber3=[ber3,ber30];
berreal=[berreal,berreal0];
end
figure
subplot(1,2,1)
k=0:5:30
semilogy(k,ber1,'-b*')
grid on
hold on
semilogy(k,ber2,'-bo')
semilogy(k,ber3,'-rx')
semilogy(k,berreal,'-.r')
axis([0 30 10^(-4) 1])
xlabel('SNR/dB'),ylabel('BER')
legend('LS','多项式','MST','理想')
set (gcf,'color',[1 1 1])
subplot(1,2,2)
k=0:5:30
semilogy(k,mse1,'-b*')
grid on
hold on
semilogy(k,mse2,'-bo')
semilogy(k,mse3,'-rx')
axis([0 30 10^(-4) 1])
xlabel('SNR/dB'),ylabel('MSE')
legend('LS','多项式','MST')
toc