www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/test_gpsArray_final_stap_hpp.m
% *** lanruru 2009.9.23
clc;
clear all;
close all;
fs=6e6;
N=fs*0.01;
tic
for i=1:6 %
fid=fopen('D:\1112界面7\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('D:\1112界面7\sv1.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1); %%% 0 改为 -1 %%%%%%
dat1y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('D:\1112界面7\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('D:\1112界面7\sv6.dat');
statusy2=fseek(fid,2*N*(i-1),-1) ; % cof or "0"%%% 0 改为 -1 %%%%%%
dat6y2=fread(fid,N,'uint16=>single'); % inf: end of the file ; 卫星3
fclose(fid)
fid=fopen('D:\1112界面7\sv6.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%%% 0 改为 -1 %%%%%%
dat6y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('D:\1112界面7\sv6.dat');
status=fseek(fid,2*N*(i-1)+4,-1);%%% 0 改为 -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));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('D:\1112界面7\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('D:\1112界面7\sv14.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%%% 0 改为 -1 %%%%%%
dat14y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('D:\1112界面7\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('D:\1112界面7\sv22.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('D:\1112界面7\sv22.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);
dat20y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('D:\1112界面7\sv22.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('D:\1112界面7\sv25.dat');
statusy2=fseek(fid,2*N*(i-1),-1); % cof or "0"
dat25y2=fread(fid,N,'uint16=>single'); % 卫星8
fclose(fid)
fid=fopen('D:\1112界面7\sv25.dat');
statusy1=fseek(fid,2*N*(i-1)+2,-1);%延迟信号
dat25y1=fread(fid,N,'uint16=>single');
fclose(fid)
fid=fopen('D:\1112界面7\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));
fc3=4.309e6;
dt=1/fs;
t=(0:N-1)*dt;
s_jam3y2=exp(1i*2*pi*fc3*t)';% sin(2*pi*fc3*t)';
guiyiJam3y2=(s_jam3y2-mean(s_jam3y2))/sqrt(var(s_jam3y2));
%%%j1的延迟
t=(1:N)*dt;
s_jam3y1=exp(1i*2*pi*fc3*t)';
guiyiJam3y1=(s_jam3y1-mean(s_jam3y1))/sqrt(var(s_jam3y1));
%%%j1的延迟
t=(2:N+1)*dt;
s_jam3=exp(1i*2*pi*fc3*t)';
guiyiJam3=(s_jam3-mean(s_jam3))/sqrt(var(s_jam3));
%色散多径
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MM=4; %天线数
TT=3; %抽头数
fs=6e6; %采样率
Ndian=fs*0.01; %采样点数 20ms长度
f0=4.309e6; %中频
fdelt=2e6; %频偏 带宽
thetaj=-50/180*pi; %直射干扰角度 -10度
thetajmp=-30/180*pi; %干扰多径角度 -50度 色散多径方向变化否？？？？？
d=0.5; %阵元间距
lambda=0.1904; %波长 L1波段
tao=1/fs;%200e-9; %时延
omiga=2*pi*f0; %中心角频率
deltomega=2*pi*fdelt; %角频率偏移
% taoi= lambda/2*sin(thetaj)/3e8; %天线间延迟时间
delaymp=150e-9; %多径时延
alpha=0.68; %多径衰减系数
p=1e-4;
D =[0:1:MM-1];
aj=exp(-1i*pi*D*sin(thetaj)).';
ajmp=exp(-1i*pi*D*sin(thetajmp)).';
%%%%%%%%%% 色散多径产生 %%%%%%%%%
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);
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
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;
RLinv=inv(RL);
Bvec=p*Rvec*RLinv;
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);
uu=disnoise(:,1:TT-1);
uu=reshape(uu,MM*(TT-1),1);
%初始化
ssss=cholsky(RL,MM*(TT-1));
disjammer0=ssss*uu;
% read_fid_real = fopen( 'e:\disjammer0.dat', 'w');
% [cou] = fwrite(read_fid_real, real(disjammer0.'), 'double');
% [cou] = fwrite(read_fid_real, imag(disjammer0.'), 'double');
% fclose(read_fid_real);
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)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%色散多径产生 完
% 构成阵列数据
% snr=-20; % dB
% jnr1=20;
% As=10^(snr/20); % 以噪声为参考，设为1
% Aj1=10^(jnr1/20);
snr=-20;
jnr1=20;
jnrm=10;
An=sqrt(p/10^(jnrm/10));
As=An*sqrt(10^(snr/10));
Aj1=An*sqrt(10^(jnr1/10));
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'];
ssj3=[Aj1*guiyiJam3';Aj1*guiyiJam3y1';Aj1*guiyiJam3y2'];
a2=-5;
a3=22;
a4=33;
a8=10;
a66=-20;
adj3=-70;
b=exp(-1i*pi*[0:MM-1]'*sin(a2*pi/180));
c=exp(-1i*pi*[0:MM-1]'*sin(a3*pi/180));
d=exp(-1i*pi*[0:MM-1]'*sin(a4*pi/180));
h=exp(-1i*pi*[0:MM-1]'*sin(a8*pi/180));
h6=exp(-1i*pi*[0:MM-1]'*sin(a66*pi/180));
doa_J3=exp(-1i*pi*[0:MM-1]'*sin(adj3*pi/180)); % 干扰3
Ab=kron(eye(3),b);
Ac=kron(eye(3),c);
Ad=kron(eye(3),d);
Ah=kron(eye(3),h);
Ah6=kron(eye(3),h6);
Aj3=kron(eye(3),doa_J3);
randn('state',100);%111 s1 101 s2
Noise=An/sqrt(2)*randn(MM*TT,N)+1i*An/sqrt(2)*randn(MM*TT,N);
x=Ab*ss1+Ac*ss14+Ad*ss20+Ah*ss25+Ah6*ss6+Noise+sjammer+Aj3*ssj3;
fid=fopen('E:\realD4_snr20_s7.dat','a');
count=fwrite(fid,real(x.'),'single');
fclose(fid)
fid=fopen('E:\imagD4_snr20_s7.dat','a');
count=fwrite(fid,imag(x.'),'single');
fclose(fid)
num
end
toc