simulation.m

利用自己仿真产生的卫星数据

tic
clear all;clc;

global datainf;

global fs
global fca;
global IF;
global num;
global datanum;

fs=16367667
ts=1/fs;
IF=20491635; 
fo=IF-1*fs; 
time=1;
datanum=fix(fs/1000);
num=[0:datanum-1];
data=signal(IF,2972,3,2760);
data1=[data data*(-1) data*(-1) data data data*(-1) data data*(-1) data*(-1) data data data*(-1) data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
data=signal(IF,1872,10,2360);
data2=[data*(-1) data data*(-1) data*(-1) data data data*(-1) data data*(-1) data*(-1) data data data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
data=signal(IF,3972,11,-4760);
data3=[data data data*(-1) data data*(-1) data data data data*(-1) data data*(-1) data data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
data=signal(IF,5172,25,4960);
data4=[data -data data -data data -data data -data data -data data -data data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
data=signal(IF,6197,28,3760);
data5=[-data data -data -data data -data -data data -data -data data -data data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
data=signal(IF,8192,30,-4960);
data6=[data data data data data data data data data data data data data data*(-1) data data*(-1) data*(-1) data data data*(-1)];
%data=data1+data2+data3+data4+data5+data6;
data=data1+data2+data3+data4+data5+data6+(rand(1,400*datanum)-0.5)*18;

Dopler=[2760 2360 -4760 4960 3760 -4960]; 
CAshift=[2972 1872 3972 5172 6197 8192];
channel=struct('satelliteID',0,...       %定义通道结构：卫星号，捕获频率，多普勒频移，码相位等信息
    'Acqfrequence',0,...
    'RealDoppler',0,...
    'Doppler1',0,...
    'Doppler2',0,...
    'Doppler3',0,...
    'AcqCAphase',0,...
    'RealPhase',0,...
    'maxAmp',0,...
    'status','0');
channel=[channel channel channel channel channel channel channel channel channel channel channel channel];
acqdata=data(1:datanum);
facqdata=fft(acqdata);
maxno=21;
maxno1=5;
k=0;
out=zeros(maxno,datanum);
for svID=3,
    fprintf('.');
    CAsample=sampleCAcode(svID,fs,time); %**按照卫星号，采样频率，和采样时间采样指定的CA码
    CAsample2=[CAsample CAsample];
    fca=fft(CAsample); 
    for fideta=1:1:21,                                  %**fideta是多普勒频移
        fi(fideta)=fo-10*1000+(fideta-1)*1000;
        carrier=exp(j*2*pi*fi(fideta)*num*ts);
        cosine= real(carrier);                                %本地信号
        sine  =imag(carrier);                                %本地信号相移90°
        I =acqdata.*cosine;   %信号也分了两相
        Q =acqdata.*sine;
        fiq=fft(I+Q*j);%此后的fft
        out(fideta,:)=abs(ifft( (fca.* conj(fiq)))).^2;      %做相关 Out是21*16367矩阵，行是频率，列代表码相位
    end
      
        [maxAmp(svID) freqshiftNo]=max(max(out'));          %每颗星的21个偏移频率中产生最大相关值的那个偏移频率值
        [maxAmp(svID) startpointNo]     =max(max(out ));    %每颗星的信号的8184 个频率点,在同21个频率做完相关后,能产生最大的相关值的那个起始位置点
        temp=out(freqshiftNo,:);
        temp(startpointNo-2:startpointNo+2)=[];
        
        
        if maxAmp(svID)/mean(temp)>37, 
            CAsample1=[CAsample(startpointNo:datanum) CAsample(1:startpointNo-1)];
            dopshift1=(freqshiftNo-1)-fix(maxno/2);
            fic=fo+dopshift1*1000;    
            stripcode=CAsample1.*acqdata;
            
            for counter=1:1:maxno1,      %-1000到1000的多普勒频移,精度为50hz
                dopshift=(counter-1)-fix(maxno1/2);
                fi=fic+dopshift*400;    %加多普勒偏移后的频率
                carrier=exp(j*2*pi*fi*num*ts);
                cosine= real(carrier);
                sine   =imag(carrier);
                I   =stripcode.*cosine;   %信号也分了两相
                Q =stripcode.*sine;
                squaresum(counter)=sum(I').^2+sum(Q').^2;
            end
            [maxAmp1 dopno]=max(squaresum);%dopno为最大的哪个频移
            k=k+1;
            channel(k).satelliteID=svID;
            channel(k).Acqfrequence=fic+((dopno-1)-fix(maxno/2))*400;
            channel(k).RealDoppler=Dopler(k);
            channel(k).RealPhase=CAshift(k);
            channel(k).Doppler1=dopshift1*1000;
            channel(k).Doppler2=((dopno-1)-fix(maxno1/2))*400+dopshift1*1000;
            channel(k).AcqCAphase=startpointNo;
            channel(k).maxAmp=maxAmp(svID);
            channel(k).maxAmpN=maxAmp1;
            
             fr=fic+((dopno-1)-fix(maxno1/2))*400;
        inputdata=data(1:8*datanum);
        CAsample10=[CAsample1 CAsample1 CAsample1 CAsample1 CAsample1 CAsample1 CAsample1 CAsample1];
        dfreq=finefrequency(CAsample10,fr,inputdata,8,k);
        accuratefreq=fr-dfreq;
        channel(k).Doppler3=channel(k).Doppler2-dfreq;
        channel(k).accuratefreq=accuratefreq;
      
        trackdata=data(9*datanum+1:10*datanum);
         %CAsample=sampleCAcode(svID,fs,time); %**按照卫星号，采样频率，和采样时间采样指定的CA码
         %CAsample2=[CAsample CAsample];
        CAsample1=CAsample2(CAshift(k):datanum+CAshift(k)-1);
        for counter=1:19,
            if counter<9,
                CAsampleP=[CAsample1(datanum-9+counter:datanum) CAsample1(1:datanum-10+counter)];
            elseif counter==9,
                CAsampleP=CAsample1;
            else
                CAsampleP=[CAsampleP(counter-9+1:datanum) CAsampleP(1:counter-9)];
            end
            CAsampleL=[CAsampleP(9:16367) CAsampleP(1:8)];
            CAsampleE=[CAsampleP(16360:16367) CAsampleP(1:16359)];
            carrier=exp(j*2*pi*fr*num*ts);
            cosine=real(carrier);
            sine  =imag(carrier);
            trackdataI=cosine.*trackdata;
            trackdataQ=sine.*trackdata;
            
           % stripcodeP=trackdata.*CAsampleP;
           % stripcodeE=trackdata.*CAsampleE;
           %  stripcodeL=trackdata.*CAsampleL;
            IP(counter)=sum(trackdataI.*CAsampleP);
            QP(counter)=sum(trackdataQ.*CAsampleP);
            IE(counter)=sum(trackdataI.*CAsampleE);
            QE(counter)=sum(trackdataQ.*CAsampleE);
            IL(counter)=sum(trackdataI.*CAsampleL);
            QL(counter)=sum(trackdataQ.*CAsampleL);
            D(counter)=((IE(counter)^2+QE(counter)^2)^0.5-(IL(counter)^2+QL(counter)^2)^0.5)/((IE(counter)^2+QE(counter)^2)^0.5+(IL(counter)^2+QL(counter)^2)^0.5);
           % D(counter)=(IE(counter)^2+QE(counter)^2)-(IL(counter)^2+QL(counter)^2);
           %D(counter)=(IE(counter)-IL(counter))*IP(counter)-(QE(counter)-QL(counter))*QP(counter);
            %n=fix(16*D(counter+1));
            %CAsample3=[CAsampleP CAsampleP CAsampleP];
            %CAsampleP=CAsample3(datanum+n:2*datanum+n-1);
           
        end
        regc(k)=min(abs(D));
        h=figure;
        subplot(4,1,1);
        plot(IP.^2+QP.^2);
        title('Present Code');
        subplot(4,1,2);
        plot(IE.^2+QE.^2);
        title('Early Code');
        subplot(4,1,3);
        plot(IL.^2+QL.^2);
        title('Late Code');
        subplot(4,1,4);
        stairs(D);
        title('D');
      % fr=fo+Dopler(k)+(-1)^k*20; 
       fr=fo+Dopler(k);
       %CAsample1=[CAsample1(3:16367) CAsample1(1:2)];
       for counter=1:395,
           trackdata=data((counter-1)*datanum+1:counter*datanum);
           carrier=exp(j*2*pi*fr*num*ts);
           cosine=real(carrier);
           sine  =imag(carrier);
           trackdataI=cosine.*trackdata;
           trackdataQ=sine.*trackdata;
           IP(counter)=sum(trackdataI.*CAsample1);
           QP(counter)=sum(trackdataQ.*CAsample1); 
       end
      h=figure;
      subplot(2,1,1);
      plot(IP);
      subplot(2,1,2);
      plot(QP);  
   end 
end
         
  
fprintf('\n');
h=figure;
stem (maxAmp, 'DisplayName', 'maxAmp', 'YDataSource', 'maxAmp');
%D=[-0.9818 -0.9874 -0.9892 -0.9930 -0.9937 -0.9942 -0.9960 -0.9963 -0.9974 -0.9831 -0.9471 -0.8755 -0.7677 -0.6048 -0.4053 -0.1530 0.1085 0.3506 0.5727 0.7409 0.8681 0.9395]