track.rar_GPS码环_GPS码跟踪环_gps锁相环_track GPS_锁相环 GPS

%实现码跟踪 % function track(carrieroff,code,codeoff) clc; clear all; format long; %%%%%%%%%%%%参数设置%%%%%%%%%%%%%%%%%%% n=4; PRNmes=[21,-500,1329;22,1500,1710;26,-3000 , 610;15,2000,101;18,500,937;6,-3500 ,536;9,3000 ,1796;3,2000,214];%PRN号，频偏，半码偏 PRN=PRNmes(n,1); carrieroff=PRNmes(n,2); codeoff=PRNmes(n,3); % 21号卫星频偏 -547 22号卫星频偏1694 26号卫星频偏-2987 9545013 15号卫星频偏1921 18号频偏246 6号卫星频偏-3688 9544312 9号卫星频偏2832 9550832 3号频偏1903 % 21号卫星偏移个半码片1329 22号卫星偏移1710个半码片 26号卫星偏移610个半码片 15号卫星偏移101个半码片 18号卫星偏移937个半码片 6号卫星偏移536个半码片 9号卫星偏移1796个半码片 3号卫星偏移214个半码片 fm=1.023e6+carrieroff/1540; fm1=1.023e6; fc=38.192e6; %采样频率 tc=1/fc; f0=1.57542*1e9; %载波频率 fi=9.548e6; %中频频率 Rb=50; %信息速率 ficar(1)=fi+carrieroff; zeta=0.707; B1=1 ; ;%0.1*Rb;%码环路噪声等效带宽 B2=9; %0.1*Rb;%载波环路噪声等效带宽 B3=40; a2=1.414; wn1=B1*4*a2/(1+a2^2);%自然角频率 1.885712990215363 wn2=B2*4*a2/(1+a2^2);%自然角频率 18.85712990215363 wn3=B3*4*a2/(1+a2^2);%自然角频率 18.85712990215363 Tc1=1e-3;%更新时间 Tc2=1e-3;%更新时间 Tc3=1e-3; %数字环路滤波器的系数 C2=wn1^2; % 3.555913481466967 C1=1.414*wn1; % 2.666398168164524 % C1=5.3429e-005; % C2=1/50*4*(wn1*Tc1)^2/(4+4*0.707*wn1*Tc1*(wn1*Tc1)^2); Cp1=wn2^3; % 6.705432244084198e+003 Cp2=1.1*wn2^2; %3.911504829613664e+002 Cp3=2.4*wn2; %45.25711176516872 Cf1=wn3^2; Cf2=1.414*wn3; codelength=2046; %%%%%%%%%%%%%%%产生捕获后的本地码相位%%%%%%%%%%%%%% CAcode1=cacode(PRN); CAcode2=CAcode1(floor(0:0.5:1022)+1); CAcode2(codelength)=CAcode1(1023); CAcode2=[CAcode2(codeoff:codelength) CAcode2(1:codeoff-1)]; if mod(codeoff,2)==1 CAcode=[CAcode2(codeoff-1) CAcode2(codeoff-1) CAcode2 CAcode2(codeoff) CAcode2(codeoff) ]; else CAcode=[CAcode2(codeoff-2) CAcode2(codeoff-1) CAcode2 CAcode2(codeoff) CAcode2(codeoff+1) ]; end %初始化变量 remCodePhase=0; Fcw_offp(1)=0; Fcw_off(1)=0; X(1)=0; Z(1)=0; Y(1)=0; df1(1)=0; df3(1)=0; pf1(1)=0; pf4(1)=0; DF(1)=0; k=1; phase(1)=0; CodeFcw_off(1)=0; while(1) codePhaseStep = 2*fm/ fc; %每个采样周期所采的码片数 blksize = ceil((fm/1000*2-remCodePhase) / codePhaseStep); fid = fopen('D:\MATLAB701\work\GPSdata-DiscreteComponents-fs38_192-if9_55.bin', 'r'); %读入接收数据 x = fread(fid,blksize*2, 'int8'); ca=x'; fclose(fid);%关闭指针 tcode=(remCodePhase-1) : ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase-1); tcode2=floor(tcode)+3; E=CAcode(tcode2); %超前即时码半个码片 tcode=(remCodePhase+1) : ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase+1); tcode2=floor(tcode)+3; L=CAcode(tcode2); %滞后即时码半个码片 tcode=remCodePhase: ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase); tcode2=floor(tcode)+3; P=CAcode(tcode2); %即时码片 remCodePhase = ((tcode(blksize) +codePhaseStep) -2046); %码相位偏移 why *0.5 +codePhaseStep trigarg=(2*pi*ficar(k)*((0:blksize*2) )*tc)+phase(k); %载波相位 dataI1=cos(trigarg(1:blksize)); %同相载波 dataQ1=sin(trigarg(1:blksize)); %正交载波 dataI2=cos(trigarg(blksize+1:blksize*2)); %同相载波 dataQ2=sin(trigarg(blksize+1:blksize*2)); %正交载波 %%%%%%%%%%%%%%%码跟踪%%%%%%%%%%%%%%%%%%% I1=dataI1.*ca((1:blksize)); I2=dataI2.*ca(blksize+1:blksize*2); Q1=dataQ1.*ca((1:blksize)); Q2=dataQ2.*ca(blksize+1:blksize*2); %%%%%%%%%%%相关积分%%%%%%%%%%%%%%%%%%%% IES(k)=sum(I1.*E); ILS(k)=sum(I1.*L); IPS1(k)=sum(I1.*P); IPS2(k)=sum(I2.*P); QES(k)=sum(Q1.*E); QLS(k)=sum(Q1.*L); QPS1(k)=sum(Q1.*P); QPS2(k)=sum(Q2.*P); % DF(k+1)=(IPS1(k)*QPS2(k)-IPS2(k)*QPS1(k)); DF(k+1)=atan2((IPS1(k)*IPS2(k)+QPS1(k)*QPS2(k)),(IPS1(k)*QPS2(k)-IPS2(k)*QPS1(k)))/(2*pi)/360*1000; IQ(k+1)=atan(QPS1(k)/IPS1(k))/(2*pi);%载波鉴相 E1(k)=sqrt(IES(k)^2+QES(k)^2); L1(k)=sqrt(ILS(k)^2+QLS(k)^2); DLL(k+1)=0.5*(E1(k)-L1(k))/(E1(k)+L1(k)); %码环鉴相误差 k=k+1; %%%%%%%%%%%%%%%%%%%%%环路滤波器鉴频%%%%%%%%%%%%%%%%%%%%%% % 二阶FLL % Z(k)=Cf2*Tc2*IQ(k)+Z(k-1); % Fcw_off1(k) =IQ(k)*Cf1+(Z(k-1)+Z(k))*0.5; Fcw_off(k)=IQ(k)-IQ(k-1); %FLL辅助的PLL pf1(k)=IQ(k)*Cp1*Tc3+Fcw_off(k-1)*Cf1*Tc2+pf1(k-1);% *4 pf2=pf1(k)+pf1(k-1); pf3=0.5*pf2+Cf2*Fcw_off(k-1)+Cp2*IQ(k); %*2 pf4(k)=pf3*Tc3+pf4(k-1); pf5=pf4(k)+pf4(k-1); P_off(k)=0.5*pf5+Cp3*IQ(k); ficar(k)=ficar(k-1)-P_off(k); phase(k)=mod(ficar(k)*2*pi/1000,2*pi)+phase(k-1); %%%%%%%%%%%%%%%%%%%%%%码环滤波器%%%%%%%%%%%%%%%%%% %码环鉴别 X(k)=C2*Tc1*DLL(k)+X(k-1); CodeFcw_off(k) =DLL(k)*C1+(X(k-1)+X(k))*0.5; % X(k)=DLL(k)*C2+X(k-1); % CodeFcw_off(k)=DLL(k)*C1+X(k); fm=fm-CodeFcw_off(k)-P_off(k)/1540 % remCodePhase=remCodePhase-CodeFcw_off(k); % fm=fm1+(ficar(k)-fi)/1540 end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot(ficar); title('ficar'); figure; plot(IPS1); title('IPS1'); figure; plot(DLL); title('DLL'); figure; plot(phase); title('phase'); %a_wave=(0:2^24:2^32-1); % memoryI=cos(2*pi/2^32*a_wave+phase);%余弦值存储rom表 % memoryQ=[memoryI(3*2^6+1:2^8),memoryI(1:3*2^6)];%移相pi/2得到正弦表 % b=mod((0:codelength-1)*(Fcw_norm-Fcw_off(k)),2^32);%将50000个点乘以频率步进后，对2^32求余 Fcw_off(k) % index=find(b>=0&b<=a_wave(1)); % dataI(index)=memoryI(1); % dataQ(index)=memoryQ(1);%找出b中从0到a_wave(1)的数据，将 dataI与dataQ中和b相应位置的数据的数据换为相应的正余弦值 % for i=2:2^8 % index=find(b>a_wave(i-1)&b<=a_wave(i)); % dataI(index)=memoryI(i); % dataQ(index)=memoryQ(i);%找出b中从a_wave(1)到a_wave(256)的数据，将 dataI与dataQ中和b相应位置的数据的数据换为相应的正余弦值 % end % index=find(b>a_wave(256)); % dataI(index)=memoryI(256); % dataQ(index)=memoryQ(256);%找出b中从大于a_wave(256)的数据，将 % dataI与dataQ中和b相应位置的数据的数据换为相应的正余弦值