GPS捕获-Matlab代码

Posted 2021-05-12 不知名的小咸鱼

1 GPS捕获原理-简要

• 首先是正交混频，剥离载波，注意本地NCO（不懂的话可以查查DDS）生成的相位步进以及相位累加器的溢出处理，没处理好会导致无法出现相关峰值，CA码不做介绍（这个看谢钢那本GPS原理和接收机设计这本书就可以了），这个没什么技术难度。然后，将CA码重采样，这里重采样为16倍码率，由于只是理论验证，且这里相关耗时并不算多，耗时最多的是正交混频和重采样。最后利用matlab的xcorr做1ms并行相关，20ms相干积分，以及寻找大于阈值的相关峰值，即可实现捕获，但由于时间较短，不能对弱信号进行捕获。
• 这里采用的是20ms的实测中频数据，参数如下： Fs = 16.386E6 Hz 中频采样率 Fc = 4.092E6
  Hz 理论中频 Fcc = 1.023E6 bps 码片速率 信号类型： GPS L1
  

Matlab代码

%% function : aquire gps signal
%% date : 20210507
%% author: skignwei
%%%%%%%%%%%%%%55555555555555555555555555%%%%%%%%%
% PRN Doppler
% 21 3740
%%%%%%%%%5555555555555555555555555%%%%%%%%%%%%5555
clear
close all

filepath = 'OutCACodeResult/';
filename1 = 'GPSCACode.txt';
%% parameter
Fs = 16.368E6;  %采样率
Fc = 4.092E6;   %理论中频数据
% Fs = 30.69E6;  %采样率
% Fc = 8.184E6;   %理论中频数据

Fcc = 1.023E6;    %码率
GPS_CACODE_NUM = 1023;
initDopp = -500;  %Hz
slew_count = Fs/Fcc; %%本地码率倍数
FFT_N = 4092;
FFT_Len = Fs/1E3;

%% read data
gpsIFdata = importdata('gps_data_20ms.mat');
gpsIFdata = gpsIFdata';
% gpsIFdata = importdata('IFdata3.txt',',');

f1 = [filepath,filename1];
d1= importdata(f1); 
%% init parameter
CarrierLoop(Fs,Fc);
Set_CarrierFreq(initDopp,Fc,Fs);
CodeLoop(19,d1);     %GPS19
Move_Chips(0);
%%
L = length(gpsIFdata);
temp_Coh_PRN = zeros(FFT_Len*2-1,32); %ALLGPS
FFTCode_Signal = zeros(FFT_Len,1);
FFTCode_CA = zeros(FFT_Len,1);
IFData = complex(gpsIFdata(:,1),0);
global reg_aquired regcount
reg_aquired = zeros(200,2);
regcount = 1;
%%并行码捕获
fcount = 1;
for m = -37.4   %%-4000Hz~4000Hz
for j = 21
% if j== 21 || j==15 || j==27 || j==9
%         continue
% end
% init parameter
tic
    FFT_count = 1;
    FFT_Cohcount = 1;
    CarrierLoop(Fs,Fc);
    initDopp = m*100;
    Set_CarrierFreq(initDopp,Fc,Fs);
    CodeLoop(j,d1);     %GPS 1-32
    Move_Chips(0);
    
    temp_Coh = zeros(FFT_Len*2-1,1);
    for i = 1:L

%     IFData = complex(gpsIFdata(i,1),gpsIFdata(i,2));

    %载波剥离,暂时不考虑多普勒
    ComplexData = ComplexDDC(IFData(i));
    %%
    %码剥离,暂时不考虑多普勒
    [tempchar,Code_Signal, Code_CA] = Code_Correlate(ComplexData,slew_count);
    if tempchar == 1
        FFTCode_Signal(FFT_count,1) = Code_Signal;
        FFTCode_CA(FFT_count,1) = complex(Code_CA,0);
    
        if FFT_count < FFT_Len
            FFT_count= FFT_count +1;
        else 
           FFT_count = 1;
    %        tic
           tempFFTCode_outVec =  xcorr(FFTCode_Signal,FFTCode_CA,'coeff'); %归一化结果，并行码相关
    %        toc 
           %%相干积分
           temp_Coh = temp_Coh + tempFFTCode_outVec;

           FFT_Cohcount = FFT_Cohcount + 1;
           if FFT_Cohcount >= 21
               disp("处理结束")
               disp(j)
           end
    %        break;
        end    
    end
    end
    temp_Coh_PRN(:,j) = temp_Coh;
toc
end
% result
temp_Coh_PRNm = abs(temp_Coh_PRN);

%%search
aquiredresult = aquire_Prn(initDopp,temp_Coh_PRNm);

%figure
figure(fcount)
fcount= fcount + 1;
plot(temp_Coh_PRNm,'DisplayName','temp_Coh_PRNm')
disp("inidopp")
disp(initDopp)
end
%% func: 载波剥离
%% return : 正交剥离后的结果

%%//-----------------------中频数据下变频----------------------------//
function Complexdata = ComplexDDC(IFData)
    global costable sintable
    global  Carrier_phase Carrier_phaseStep
	%%Multiply 小心这里没处理好
	Carrier_phase = Carrier_phase + Carrier_phaseStep;			%相位步进
    if Carrier_phase >= pow2(32)
       Carrier_phase = mod(Carrier_phase,pow2(32));
    end
	phaseIndex = bitshift(Carrier_phase,-27,'uint32');
    phaseIndex = mod(phaseIndex,32)+1;
    %index is different to cpp
	TempComplexData = complex(costable(phaseIndex), sintable(phaseIndex));						%%正频率
	TempComplexData = IFData * TempComplexData;
	TempComplexData = complex(bitshift(real(TempComplexData),-4,'int64'), bitshift(imag(TempComplexData),-4,'int64'));
   
    %%return
    Complexdata = TempComplexData;


end

%% //-----------------------载波频率设定----------------------------//
function Set_CarrierFreq(DopplerFreq,GPS_FREQ_COMPENSATE,SystemClockFre)
    global  Carrier_phaseStep
	double tempData;
	tempData = DopplerFreq + GPS_FREQ_COMPENSATE;
	tempData = pow2(32) * tempData / SystemClockFre;
	Carrier_phaseStep = floor(tempData);
end

%% //-----------------------初始 载波频率设定----------------------------//
function CarrierLoop(SystemClockFre,GPS_FREQ_COMPENSATE)
    global costable sintable
    global  Carrier_phase Carrier_phaseStep
    %% 正弦波表
    costable = [252, 247, 233, 210, 178, 140, 96, 49, 0, -49, -96, -140, -178, -210, -233, -247, -252, -247, -233, -210, -178, -140, -96, -49, 0, 49, 96, 140, 178, 210, 233, 247 ];
    sintable = [ 0, 49, 96, 140, 178, 210, 233, 247, 252, 247, 233, 210, 178, 140, 96, 49, 0, -49, -96, -140, -178, -210, -233, -247, -252, -247, -233, -210, -178, -140, -96, -49 ];

	double tempData;
	Carrier_phase = 0;
	%% 计算标称中频所对应的M
	tempData = GPS_FREQ_COMPENSATE;
	tempData = pow2(32) * tempData / SystemClockFre;
	Carrier_phaseStep = round(tempData);
    
end

%% func: 码剥离
%% return : 码剥离后的结果
%% 并行相关器
%% result 相关结果
function [temp_char,FFTCode_Signal, FFTCode_CA]= Code_Correlate(DDCData,slew_count)
    global Code_phaseIndex
    global CA_code count_cc 
    global div_count
    div_count = div_count + 1;
    if div_count>0
        div_count =0;
        FFTCode_Signal = DDCData;
        %%重采样
        FFTCode_CA = CA_code(Code_phaseIndex);
        if count_cc < slew_count
             count_cc =  count_cc + 1;
        else
             count_cc = 1;
             Code_phaseIndex = Code_phaseIndex+1;        
             if Code_phaseIndex >1023
                 Code_phaseIndex = 1;
             end

        end
        temp_char = 1;
    else
        temp_char = 0;   
        FFTCode_Signal = 0;
        FFTCode_CA = 0;
    end

end


%% 码片搜索设定
function Move_Chips(chip)
    global CurrentCodePhase slewCount
	tempdata = chip - CurrentCodePhase;
	CurrentCodePhase = chip;

    if tempdata >= 0
        slewCount = tempdata;
    else
        slewCount = 16 * GPS_CACODE_NUM + tempdata;   %16倍码率
    end
end

%% 初始设定
function CodeLoop(PRN,CA)
    global CurrentCodePhase   Code_phaseIndex 
    global CA_code count_cc
    global div_count
    div_count =0;
    CurrentCodePhase = 0;
    Code_phaseIndex = 1;
    count_cc = 1;
    CA_code = CA(:,PRN);
end

%% 捕获结果
function aquiredresult = aquire_Prn(dopp,recorr)
 global reg_aquired regcount
 b = max(recorr,[],1);  %%按列
 index = find(b>0.3);  %%捕获阈值
 if isempty(index)
 else
    L =length(index);
    PRN = index;
    initdopp = dopp;
    for i = 1:L
        reg_aquired(regcount,:) = [PRN(i),initdopp];
        regcount = regcount+1;
    end   
 end
 aquiredresult = reg_aquired;
end

捕获结果

这里给出GPS21的捕获结果，因为matlab的for循环是很耗时间的，为了逻辑上和实际近似，没有采用矩阵运算，所以做一次捕获费时12s（我的电脑配置还算好的）。后面就懒得移植C++了，毕竟我要画图，画图，C++用多线程估计1秒不到就可以处理晚了吧。我试了一下，大概4小时完成一次GPS1-32的并行捕获，多普勒范围是-4000-4000Hz，步进200Hz，其实100Hz会好一些。

数据

想要数据可以私戳我，我个人不太喜欢传到网盘