www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/test_gpsArray_final_stap_twomp_three.m
% ***** 由单颗卫星数据合成单天线数据（>=4颗），然后生成阵列天线数据，
% *** lanruru 2009.9.23
clc;
clear all; close all;
% 读取数据
fs=5e6;%20e6;
N=fs*0.02; %
p=1e-6; % p 为色散多径干扰功率
snr=-20;
jnr1=25; %点频干扰干躁比
jnrm=25; %色散多径干扰干躁比
An=sqrt(p/10^(jnrm/10));
As=An*sqrt(10^(snr/10));
Aj1=An*sqrt(10^(jnr1/10));
tic
for i=1:1850 % 1 s data length，3690 3700
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv1.dat');
statusy2=fseek(fid,2*N*(i-1),-1); % cof or "0" %%% 0 改为 -1 %%%%%%
dat1y2=fread(fid,N,'uint16=>single'); % inf: end of the file ; 卫星2
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv1.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1); %%% 0 改为 -1 %%%%%%
dat1y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv1.dat');
status=fseek(fid,2*N*(i-1)+4,-1);%%% 0 改为 -1 %%%%%%
dat1=fread(fid,N,'uint16=>single');
fclose(fid)
guiyiDat1=(dat1-mean(dat1))/sqrt(var(dat1));
guiyiDat1y1= (dat1y1-mean(dat1y1))/sqrt(var(dat1y1));
guiyiDat1y2= (dat1y2-mean(dat1y2))/sqrt(var(dat1y2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv14.dat');
statusy2=fseek(fid,2*N*(i-1),-1) ; % cof or "0"%%% 0 改为 -1 %%%%%%
dat14y2=fread(fid,N,'uint16=>single'); % inf: end of the file ; 卫星3
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv14.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%%% 0 改为 -1 %%%%%%
dat14y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv14.dat');
status=fseek(fid,2*N*(i-1)+4,-1);%%% 0 改为 -1 %%%%%%
dat14=fread(fid,N,'uint16=>single');
fclose(fid)
guiyiDat14=(dat14-mean(dat14))/sqrt(var(dat14));
guiyiDat14y1= (dat14y1-mean(dat14y1))/sqrt(var(dat14y1));
guiyiDat14y2= (dat14y2-mean(dat14y2))/sqrt(var(dat14y2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv20.dat');
statusy2=fseek(fid,2*N*(i-1),-1) % cof or "0"
dat20y2=fread(fid,N,'uint16=>single'); % inf: end of the file ; 卫星4
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv20.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);
dat20y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv20.dat');
status=fseek(fid,2*N*(i-1)+4,-1);
dat20=fread(fid,N,'uint16=>single');
fclose(fid)
guiyiDat20=(dat20-mean(dat20))/sqrt(var(dat20));
guiyiDat20y1= (dat20y1-mean(dat20y1))/sqrt(var(dat20y1));
guiyiDat20y2= (dat20y2-mean(dat20y2))/sqrt(var(dat20y2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv25.dat');
statusy2=fseek(fid,2*N*(i-1),-1); % cof or "0"
dat25y2=fread(fid,N,'uint16=>single'); % 卫星8
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv25.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%延迟信号
dat25y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv25.dat');
status=fseek(fid,2*N*(i-1)+4,-1);%延迟信号
dat25=fread(fid,N,'uint16=>single');
fclose(fid)
guiyiDat25=(dat25-mean(dat25))/sqrt(var(dat25));
guiyiDat25y1= (dat25y1-mean(dat25y1))/sqrt(var(dat25y1));
guiyiDat25y2= (dat25y2-mean(dat25y2))/sqrt(var(dat25y2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv6.dat');
statusy2=fseek(fid,2*N*(i-1),-1); % cof or "0"
dat6y2=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv6.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%延迟信号
dat6y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\tag模拟器软件修正版\sv6.dat');
status=fseek(fid,2*N*(i-1)+4,-1);%延迟信号
dat6=fread(fid,N,'uint16=>single');
fclose(fid)
guiyiDat6=(dat6-mean(dat6))/sqrt(var(dat6));
guiyiDat6y1= (dat6y1-mean(dat6y1))/sqrt(var(dat6y1));
guiyiDat6y2= (dat6y2-mean(dat6y2))/sqrt(var(dat6y2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fc=0.75e6;%4.309e6;
dt=1/fs;%50ns
t=(0:N-1)*dt;
s_jamy2=exp(1i*2*pi*fc*t)';%sin(2*pi*fc*t)';
guiyiJamy2=(s_jamy2-mean(s_jamy2))/sqrt(var(s_jamy2));
%%%j1的延迟
t=(1:N)*dt;
s_jamy1=exp(1i*2*pi*fc*t)';%sin(2*pi*fc*t)';
guiyiJamy1=(s_jamy1-mean(s_jamy1))/sqrt(var(s_jamy1));
%%%j1的延迟
t=(2:N+1)*dt;
s_jam=exp(1i*2*pi*fc*t)';%sin(2*pi*fc*t)';
guiyiJam=(s_jam-mean(s_jam))/sqrt(var(s_jam));
ssj2=[Aj1*guiyiJam';Aj1*guiyiJamy1';Aj1*guiyiJamy2'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MM=4; %天线数
TT=3; %抽头数
fs=5e6; %采样率
Ndian=fs*0.02; %采样点数 20ms长度
f0=0.55e6; %中频
fdelt=2e6; %频偏 带宽
thetaj=-70/180*pi; %直射干扰角度
thetajmp=85/180*pi; %干扰多径角度
thetajmp2=70/180*pi;
d=0.5; %阵元间距
lambda=0.1904; %波长 L1波段
tao=200e-9; %%100*10^-9; %时延
omiga=2*pi*f0; %中心角频率
deltomega=2*pi*fdelt; %角频率偏移
% taoi= lambda/2*sin(thetaj)/3e8; %天线间延迟时间
delaymp=150e-9; %60ns的多径时延
delaymp2=400e-9;
delaymp_=delaymp2-delaymp;
alpha=0.88; %多径衰减系数
alpha2=0.68;
D =[0:1:MM-1];
aj=exp(-1i*pi*D*sin(thetaj)).';
ajmp=exp(-1i*pi*D*sin(thetajmp)).';
ajmp2=exp(-1i*pi*D*sin(thetajmp2)).';
%%%%%%%%%%%%%%% 色散多径产生 %%%%%%%%%%%%%
taoT=[0 tao 2*tao;
-tao 0 tao;
-2*tao -tao 0];
rT=sinc(deltomega.*taoT./2).*exp(1i.*omiga.*taoT);
taot=[delaymp delaymp+tao delaymp+2*tao;
delaymp-tao delaymp delaymp+tao;
delaymp-2*tao delaymp-tao delaymp];
rmp=sinc(deltomega.*taot./2).*exp(1i.*omiga.*taot);
taot2=[delaymp delaymp-tao delaymp-2*tao;
delaymp+tao delaymp delaymp-tao;
delaymp+2*tao delaymp+tao delaymp];
rmp2=sinc(deltomega.*taot2./2).*exp(1i.*omiga.*taot2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
taot_m2=[delaymp2 delaymp2+tao delaymp2+2*tao;
delaymp2-tao delaymp2 delaymp2+tao;
delaymp2-2*tao delaymp2-tao delaymp2];
rmp_m2=sinc(deltomega.*taot_m2./2).*exp(1i.*omiga.*taot_m2);
taot2_m2=[delaymp2 delaymp2-tao delaymp2-2*tao;
delaymp2+tao delaymp2 delaymp2-tao;
delaymp2+2*tao delaymp2+tao delaymp2];
rmp2_m2=sinc(deltomega.*taot2_m2./2).*exp(1i.*omiga.*taot2_m2);
%%%%%%%%%%%%%%%%%%%%%%%%
taot_m_=[delaymp_ delaymp_+tao delaymp_+2*tao;
delaymp_-tao delaymp_ delaymp_+tao;
delaymp_-2*tao delaymp_-tao delaymp_];
rmp_m_=sinc(deltomega.*taot_m_./2).*exp(1i.*omiga.*taot_m_);
taot2_m_=[delaymp_ delaymp_-tao delaymp_-2*tao;
delaymp_+tao delaymp_ delaymp_-tao;
delaymp_+2*tao delaymp_+tao delaymp_ ];
rmp2_m_=sinc(deltomega.*taot2_m_./2).*exp(1i.*omiga.*taot2_m_);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for i=1:TT
% for j=1:TT
% rtao(:,:,i,j)=p*[aj ajmp]*[rT(i,j) alpha*conj(rmp2(i,j));conj(alpha)*rmp(i,j) abs(alpha^2)*rT(i,j)]*[aj ajmp]';%归一化后再乘p
% end
% end
for i=1:TT
for j=1:TT
rtao(:,:,i,j)=p*[aj ajmp ajmp2]*[rT(i,j) alpha*conj(rmp2(i,j)) alpha2*conj(rmp2_m2(i,j));conj(alpha)*rmp(i,j) abs(alpha^2)*rT(i,j) conj(alpha)*alpha2*conj(rmp2_m_(i,j));conj(alpha2)*rmp_m2(i,j) conj(alpha2)*alpha*rmp_m_(i,j) abs(alpha2^2)*rT(i,j)]*[aj ajmp ajmp2]';%归一化后再乘p
end
end
R_TT=[rtao(:,:,1,1) rtao(:,:,1,2) rtao(:,:,1,3);
rtao(:,:,2,1) rtao(:,:,2,2) rtao(:,:,2,3);
rtao(:,:,3,1) rtao(:,:,3,2) rtao(:,:,3,3)];
R0=rtao(:,:,1,1);
RL=[rtao(:,:,1,1) rtao(:,:,1,2);
rtao(:,:,2,1) rtao(:,:,2,2)];
Rvec=[rtao(:,:,1,2) rtao(:,:,1,3)];
Bvec=p*Rvec/RL;
% C=(R0-Bvec*Rvec')/p;
Wvec=[eye(MM) -Bvec(:,1:MM) -Bvec(:,MM+1:2*MM)];%%%%%%%%随天线数变
C=Wvec*R_TT*Wvec';
A0=cholsky(C,MM);
disjammer=zeros(MM,Ndian);
randn('state',0);
disnoise=randn(MM,Ndian+TT-1);%sqrt(p)*
uu=disnoise(:,1:TT-1);
uu=reshape(uu,MM*(TT-1),1);
%初始化
disjammer0=cholsky(RL,MM*(TT-1))*uu;
for kk=1:TT-1
disjammer(:,kk)=disjammer0((kk-1)*MM+1:(kk-1)*MM+MM);
end
%迭代产生
for i=TT:Ndian+TT-1
disjammer(:,i)=Bvec*[disjammer(:,i-1);disjammer(:,i-2)]+A0*disnoise(:,i);
end
sjammer=[disjammer(:,3:Ndian+2);disjammer(:,2:Ndian+1);disjammer(:,1:Ndian)];
ss1=[As*guiyiDat1';As*guiyiDat1y1';As*guiyiDat1y2'];
ss14=[As*guiyiDat14';As*guiyiDat14y1';As*guiyiDat14y2'];
ss20=[As*guiyiDat20';As*guiyiDat20y1';As*guiyiDat20y2'];
ss25=[As*guiyiDat25';As*guiyiDat25y1';As*guiyiDat25y2'];
ss6=[As*guiyiDat6';As*guiyiDat6y1';As*guiyiDat6y2'];
a1=-25;
a2=-5;
a3=5;
a4=15;
a6=30;
a102=-60;
aa=exp(-1i*pi*[0:MM-1]'*sin(a1*pi/180));
bb=exp(-1i*pi*[0:MM-1]'*sin(a2*pi/180));
cc=exp(-1i*pi*[0:MM-1]'*sin(a3*pi/180));
dd=exp(-1i*pi*[0:MM-1]'*sin(a4*pi/180));
ff=exp(-1i*pi*[0:MM-1]'*sin(a6*pi/180));
doa_J2=exp(-1i*pi*[0:MM-1]'*sin(a102*pi/180));
Aa=kron(eye(3),aa);
Ab=kron(eye(3),bb);
Ac=kron(eye(3),cc);
Ad=kron(eye(3),dd);
Af=kron(eye(3),ff);
Aj2=kron(eye(TT),doa_J2);
randn('state',134);
Noise=An/sqrt(2)*randn(MM*TT,N)+1i*An/sqrt(2)*randn(MM*TT,N);
x=Aa*ss1+Ab*ss14+Ac*ss20+Ad*ss25+Af*ss6+Noise+sjammer+Aj2*ssj2;
% figure;set(gcf, 'color', 'white');
% plot(abs(fft(x(1,1:5000))))
%
% % Rxj=x*x'/400;
% Rxj=x(:,1:40000)*x(:,1:40000)'/400;
% rank2= sort((eig(Rxj)), 'descend');
% figure;set(gcf, 'color', 'white');
% plot(1:MM*TT,10*log10(rank2),'Linewidth',1.0);
% xlabel('特征值分布');
fid=fopen('E:\BJ一朴GPS资料\realD4_snr20.dat','a'); % 列优先存!!!w改成a
count=fwrite(fid,x.','single');
fclose(fid)
fid=fopen('E:\BJ一朴GPS资料\imagD4_snr20.dat','a');
count=fwrite(fid,imag(x.'),'single'); % 虚部 w改成a
fclose(fid)
end
toc