1.zip_channel estimation_indoor_可见光_可见光通信_可见光通信信道模型

addpath('G:\host\program\VLC MIMO-OFDM\tensor_toolbox'); clc;clear tic for z=1:20 Ns=4;%LED数 Nr=4;%PD数 Nd=4;%PD数 subcarrier_num=64;%子载波数 symbol_num=32;%OFDM符号数 Ng=subcarrier_num/8; Np0=symbol_num/8;% Numbers of 0 pilots Np1=subcarrier_num/2/4*symbol_num/2-Np0; Nbps=4; M=2^Nbps; % Number of bits per (modulated) symbol %Es=7; A=sqrt(3/2/(M-1)*Es); % Signal energy and QAM normalization factor Xp1 =ones(1,Np1) ;% Pilot sequence generation data=randi([0,1],1,(subcarrier_num*symbol_num/2-symbol_num-Np1));%生成一控制命令符号序列 Data = qammod(data,M);%16qam调制 %----------------------------------插入导频------------------------------------- Xf=zeros(subcarrier_num*symbol_num/2,1);Xf(1)=1; Ps1_3=[4,4,4,4,4,4,4,37];%导频间隔 Ps4=[4,4,4,4,4,4,4,33]; for i=0:15 for j=1:8 if mod((i+1),4)==0 Psn(j+8*i)=Ps4(j); else Psn(j+8*i)=Ps1_3(j); end end end Pn=zeros(1,128); Pn(1)=1; for j=1:127 Pn(j+1)=Psn(j)+Pn(j); end Xf(Pn)=1; %插入1位置导频 for i=1:subcarrier_num*symbol_num/2 if Xf(i)==0 Xf(i)=2; end end k0=1; for p=1:32 for q=1:32 xx(q,p)=Xf(k0); Xfn(q,p)=Xf(k0);% 串并变换 k0=1+k0; end end Xfn(1,:)=0; k1=1; for p=1:32 for q=1:32 if Xfn(q,p)==2 Xfn(q,p)=Data(k1); k1=k1+1; % else % Xfn(q,p)=Xfn(q,p); end end end Xfc=flipud(conj(Xfn(2:subcarrier_num/2,:)));%子载波分配 X0=zeros(symbol_num,1); x=[Xfn' X0 Xfc']'; %导频位置定位 xxp1=x; xn=reshape(x,subcarrier_num*symbol_num,1); xp1=find(xn==1); xp0=find(xn==0); for i=1:subcarrier_num for j=1:symbol_num if xxp1(i,j)~=1 xxp1(i,j)=0; end end end %产生OFDM frame for sn=1:symbol_num X_ifft(:,sn)=ifft(x(:,sn),subcarrier_num); %ifft end Xcp =[X_ifft(subcarrier_num-Ng+1:subcarrier_num,:) ;X_ifft]; % Insert CP Xt=[Xcp Xcp Xcp Xcp]; X=reshape(Xt,(subcarrier_num+Ng),symbol_num,Ns); %P/S for t=1:subcarrier_num+Ng for tt=1:symbol_num %每路子载波上有四个导频 for ttt=1:Ns X_ALS(tt,ttt,t)=X(t,tt,ttt); end end end %------------------------------channel S-R----------------------------------- ARX = 0.01; psic = 70; %接收机视角 phi_half = 100;%半功率角度 g = -log(2)./log(cos(phi_half));%朗伯指数 Ro = (g+1)/(2*pi*(cos(pi/3))^g); for n=1:Ns for m=1:Nr if n==m psi(n,m)=atan(sqrt(2)); d(n,m)=1/cos(psi(n,m)); hSR(n,m) = ARX/(d(n,m)^2)*Ro*cos(psi(n,m)); elseif abs(m-n)==2 psi(n,m)=atan(1.5*sqrt(2)); d(n,m)=1/cos(psi(n,m)); hSR(n,m) = ARX/(d(n,m)^2)*Ro*cos(psi(n,m)); else psi(n,m)=atan(sqrt(3.25)); d(n,m)=1/cos(psi(n,m)); hSR(n,m) = ARX/(d(n,m)^2)*Ro*cos(psi(n,m)); % psi<psic end end end % %--------------------------------H_RD--------------------------------------- for a=1:Nr for b=1:Nd switch a case 1 d1=[1.5*sqrt(2),2/sin(atan(2/1.5)),2/sin(atan(2)),2/sin(atan(1.5))]; psi0(a,b)=atan(d1(b)/2.15); drd(a,b)=2.15/cos(psi0(a,b)); hRD(a,b) = ARX/(drd(a,b)^2)*Ro*cos(psi0(a,b)); case 2 d2=[1/sin(atan(1/1.5)),1.5*sqrt(2),1/sin(atan(1/1.5)),sqrt(2)]; psi0(a,b)=atan(d2(b)/2.15); drd(a,b)=2.15/cos(psi0(a,b)); hRD(a,b) = ARX/(drd(a,b)^2)*Ro*cos(psi0(a,b)); case 3 d3=[2/sin(atan(2)),2/sin(atan(2/1.5)),1.5*sqrt(2),2/sin(atan(1.5))]; psi0(a,b)=atan(d3(b)/2.15); drd(a,b)=2.15/cos(psi0(a,b)); hRD(a,b) = ARX/(drd(a,b)^2)*Ro*cos(psi0(a,b)); case 4 d4=[2/sin(atan(2/1.5)),2*sqrt(2),2/sin(atan(2/1.5)),1.5*sqrt(2)]; psi0(a,b)=atan(d4(b)/2.15); drd(a,b)=2.15/cos(psi0(a,b)); hRD(a,b) = ARX/(drd(a,b)^2)*Ro*cos(psi0(a,b)); end end end % %---------------------------------relay----------------------------------------- % relay=inv(hSR)*inv(hRD); %中继放大转发器 固定增益 % % [V Relay]=eig(relay); % %--------------------------------recieve---------------------------------------- SNR = [0:2:20]; for p=1:11 ACG=sum(sum(abs(hSR)))/16; Beta(p)=sqrt(1/(ACG^2+1./(10.^(SNR(p)/10)))); Relay{p}=Beta(p)*eye(4,4);%中继增益 % Relay=sqrt(1/(abs(hSR)^2*1+1/(10.^(SNR(p)/10))));%中继增益 for sn=1:subcarrier_num+Ng R{p,sn}=awgn(hSR'*X_ALS(:,:,sn)',SNR(p),'measured'); Ycp{p,sn}=awgn(hRD'*Relay{p}*R{p,sn},SNR(p),'measured');%S/P for tt=1:symbol_num for ttt=1:Nd ycp{p}(sn,tt,ttt)=Ycp{p,sn}(ttt,tt); end end end Y{p}=ycp{p}(Ng+1:subcarrier_num+Ng,:,:); % Remove CP for syn=1:symbol_num for ttt=1:Nd Y_fft{p}(:,syn,ttt)=fft(Y{p}(:,syn,ttt)); %%FFT end end % --------------------------------信道估计-------------------------------------------- lk=1; for nd=1:Nd for i= 1:subcarrier_num for j=1:symbol_num if xxp1(i,j)==1 K{p}(i,j,nd)=Y_fft{p}(i,j,nd); L{p}(i,lk,nd)=Y_fft{p}(i,j,nd); lk=lk+1;%1位置导频提取 if lk>4 lk=1; end end end end end L{p}(33,:,:)=[]; L{p}(1,:,:)=[]; L0{p}(1,:,:)=Y_fft{p}(1,:,:); %0位置导频提取 L0{p}(2,:,:)=Y_fft{p}(33,:,:); Yp0{p}=(L0{p}(1,:,:)+L0{p}(2,:,:))/2; for t=1:subcarrier_num-2 for tt=1:4 %每路子载波上有四个导频 for ttt=1:Nd ALS{p}(ttt,tt,t)=L{p}(t,tt,ttt); end end end [ykh,ykl]=find(xxp1'==1); klie=reshape(ykh',4,subcarrier_num-2); for t=1:subcarrier_num-2 y0weizhi=klie(:,t); for ttt=1:Nd for tt=1:4 ye{p}(ttt,tt,t)=Yp0{p}(1,y0weizhi(tt),ttt); end end end %——————————————————ELS—————————————————————— xx1=ones(Ns,4); teye=[];EstiY0{1,p}=[];Drelay{p}=[]; B{p}=[];EstY_new{p}=[];A{p}=[]; for t=1:subcarrier_num-2 teye=[teye eye(4,4)]; Drelay{p}=[Drelay{p} Relay{p}]; vecX{p,t}=reshape(ALS{p}(:,:,t)',Nd*4,1); vecye{p,t}=reshape(ye{p}(:,:,t)',Nd*4,1); end Teye{p}=reshape(teye,Nd,4,subcarrier_num-2); Trelay{p}=reshape(Drelay{p},Nd,4,subcarrier_num-2); for t=1:subcarrier_num-2 Y0{p}(:,:,t)=hRD'*Trelay{p}(:,:,t)*hSR'*Teye{p}(:,:,t)*xx1'; EstiY0{1,p}=[EstiY0{1,p} Y0{p}(:,:,t)]; end flag=0;c=1;hsr{c,p}=hSR'; while flag==0 EstiY0{c+1,p}=[];F1{p}=[];Vex{p}=[];Z2{p}=[];F4{p}=[];A{p}=[];B{p}=[];Vexe{p}=[];F3{p}=[]; Ye{p}=[];yb{p}=[]; for t=1:subcarrier_num-2 Ye{p}=[Ye{p} ye{p}(:,:,t)]; yb{p}=[yb{p};ye{p}(:,:,t)]; A{p}=[A{p} ALS{p}(:,:,t)]; vecX{p,t}=reshape(ALS{p}(:,:,t)',Nd*4,1); Vex{p}=[Vex{p};vecX{p,t}]; vecye{p,t}=reshape(ye{p}(:,:,t)',Nd*4,1); Vexe{p}=[Vexe{p};vecye{p,t}]; f1{p}(:,:,t)=Trelay{p}(:,:,t)*hsr{c,p}*Teye{p}(:,:,t)*xx1'; F1{p}=[F1{p} f1{p}(:,:,t)]; end Hrd{p}=(A{p}-Ye{p})*pinv(F1{p}); %估计rd for t=1:subcarrier_num-2 z2{p,t}=kron(xx1*Teye{p}(:,:,t),Hrd{p}*Tre