代码和测试数据: gitee link
1 GPS捕获原理
首先是正交混频，剥离载波，CA码不做介绍（参考谢钢那本GPS原理和接收机设计）。将CA码重采样，对于GPS L1CA，考虑计算量，采用2倍采样率，即2.046MHz。利用matlab的xcorr做1ms并行相关，N次相干积分，M次非相干积分，寻找大于阈值的非相干结果的峰值，即可实现捕获。
2 Matlab代码
gnss_main.m
%%
clear
close all
%% gps L1 bds B1I prn and doppler
% 4 7 8 9 16 21 26 27
% -2272 3143 777 -385 -1963 3126 -2332 -1556
% 1 2 3 4 5 6 7 9 10 16 27 30
% 10 -4 24 25 4 -1021 -138 -446 -60 -1019 -1255 1757
%% datapath and read data
ifdata_path = 'IFData\\';
codepath = 'OutCACodeResult\\';
pathname = 'ifdata_2bit_20ms.txt';
dataformat = 2; %% 1: txt IQ, with header, 2: txt no header
system_type = 1; %% 0: GPS L1CA 1: BDS: B1I
%% read raw ifdata
switch dataformat
case 1
filename = [ifdata_path,pathname];
readdata = importdata(filename);
gnss_ifdata = readdata.data;
case 2
filename = [ifdata_path,pathname];
readdata = importdata(filename);
gnss_ifdata = readdata;
end
%% read code data
switch system_type
case 0
filename1 = 'GPSCACode.txt';
f1 = [codepath,filename1];
gnss_code = importdata(f1);
case 1
filename1 = 'BDSCACode.txt';
f1 = [codepath,filename1];
gnss_code = importdata(f1);
case 2
disp("error system!!!")
case 3
disp("error system!!!")
otherwise
disp("error system!!!")
end
%% gnss config init
gnss_config = set_gnss_config(system_type);
%% acq
tic
acq_result = gnss_acq(gnss_config.acq_config, gnss_ifdata, gnss_code);
toc
if acq_result.acq_status == 1
%% show result
dop_list = gnss_config.acq_config.acq_fre_min : gnss_config.acq_config.acq_fre_step : gnss_config.acq_config.acq_fre_max;
code_phase_step = gnss_config.acq_config.code_phase_step;
system_name = gnss_config.acq_config.system_name;
x = acq_result.acq_check_list{1,1}.x * code_phase_step;
acq_noncoh_result = acq_result.acq_check_list{1,1}.temp_noncoh_result;
PRN = acq_result.acq_check_list{1,1}.prn;
%% 2D plot
figure(1)
plot(x, acq_noncoh_result(:, acq_result.acq_check_list{1, 1}.acq_check_result.dop_cnt))
xlabel('Codephase/chip')
ylabel('Amplitude')
titlename = [system_name, ' PRN-', num2str(PRN),' ','Result'];
title(titlename)
grid on
% 3D plot, note: gnss_acq not save result, please config
if 1
h = figure(2);
mesh(dop_list, x, acq_noncoh_result)
titlename = [system_name, ' PRN-', num2str(PRN),' ','Result'];
title(titlename)
xlabel('Doppler/Hz')
ylabel('Codephase/chip')
zlabel('Amplitude')
% savefig(h, titlename);
% close(h)
end
end
set_gnss_config.m
%% function: set_gnss_config
%% system_type: 0: L1CA, 1: B1I, 2: E1 3: F1
function gnss_config = set_gnss_config(system_type)
switch system_type
case 0
acq_config.system_type = 0;
acq_config.system_name = 'GPS L1CA';
acq_config.ifdata_center_fre = 1567.236E6; %Hz
acq_config.samplerate = 30.69E6;
acq_config.system_fre = 1575.42E6;
acq_config.acq_samplerate = 2.046E6;
acq_config.code_rate = 1.023E6;
acq_config.code_len = 1023;
acq_config.isdownsample = 1;
acq_config.coh_num = 1;
acq_config.noncoh_num = 20;
acq_config.acq_snr_thd = 3; % 12: 10mscoh 3: 1mscoh
acq_config.acq_prn_list = [4 7 8 9 16 21 26 27];
acq_config.is_set_acq_prn = 1;
acq_config.acq_fre_step = 1000/acq_config.coh_num/2;
acq_config.acq_fre_max = 4000;
acq_config.acq_fre_min = -4000;
acq_config.code_phase_step = 0.5;
case 1
acq_config.system_type = 1;
acq_config.system_name = 'BDS B1I';
acq_config.ifdata_center_fre = 1567.236E6; %Hz
acq_config.samplerate = 30.69E6;
acq_config.system_fre = 1561.098E6;
acq_config.acq_samplerate = 4.092E6;
acq_config.code_rate = 2.046E6;
acq_config.code_len = 2046;
acq_config.isdownsample = 1;
acq_config.coh_num = 1;
acq_config.noncoh_num = 20;
acq_config.acq_snr_thd = 3;
acq_config.acq_prn_list = [1 2 3 4 5 6 7 9 10 16 27 30];
acq_config.is_set_acq_prn = 1;
acq_config.acq_fre_step = 1000/acq_config.coh_num/2;
acq_config.acq_fre_max = 4000;
acq_config.acq_fre_min = -4000;
acq_config.code_phase_step = 0.5;
case 2
disp("no support system!!!")
case 3
disp("no support system!!!")
otherwise
disp("error system!!!")
end
gnss_config.acq_config = acq_config;
end
gnss_acq.m
%%
%% function: fft_acq for GPS L1CA, BDS B1I
%% input: gnss_config.acq_config gnss_ifdata gnss_code
%% system_type: 0: L1CA, 1: B1I, 2: E1 3: F1
%% gnss_ifdata: (real, imag) Nx2
%% gnss_code: CA_Code
%% output: acq result of all/set satellites
function acq_result = gnss_acq(acq_config, gnss_ifdata, gnss_code)
%% get config
system_type = acq_config.system_type;
ifdata_center_fre = acq_config.ifdata_center_fre;
samplerate = acq_config.samplerate;
system_fre = acq_config.system_fre;
acq_samplerate = acq_config.acq_samplerate;
code_rate = acq_config.code_rate;
code_len = acq_config.code_len;
isdownsample = acq_config.isdownsample;
%%acq parameter
coh_num = acq_config.coh_num;
acq_prn_list = acq_config.acq_prn_list;
noncoh_num = acq_config.noncoh_num;
acq_snr_thd = acq_config.acq_snr_thd;
is_set_acq_prn = acq_config.is_set_acq_prn;
acq_fre_step = acq_config.acq_fre_step;
acq_fre_max = acq_config.acq_fre_max;
acq_fre_min = acq_config.acq_fre_min;
code_phase_step = acq_config.code_phase_step;
acq_ms = coh_num*noncoh_num;
%% start deal
acq_result.acq_status = 0;
ifdata_ms = length(gnss_ifdata(:,1));
if acq_ms > ifdata_ms
acq_result.acq_status = 0;
else
ifdata_len = acq_ms/1000*samplerate;
acq_ifdata = gnss_ifdata(1:ifdata_len,:);
switch system_type
case 0
disp("ACQ GPS L1CA!!!")
%%
if is_set_acq_prn == 0
acq_prn_list = 1:32;
end
case 1
disp("ACQ BDS B1I!!!")
if is_set_acq_prn == 0
acq_prn_list = 1:37;
end
case 2
disp("no support system!!!")
case 3
disp("no support system!!!")
otherwise
disp("error system!!!")
end
% isdownsample check
if isdownsample == 1
downfre_config.fc = system_fre - ifdata_center_fre;
downfre_config.fc_down = system_fre - system_fre;
downfre_config.fs_down = acq_samplerate;
downfre_config.fs = samplerate;
% down convert fre
conv_acq_data = downconvertsample(acq_ifdata, downfre_config);
down_fre = 0;
acq_ifdata_ms_len = acq_samplerate/1000;
acq_ifdata_len = acq_ms/1000*acq_samplerate;
else
acq_samplerate = samplerate;
acq_ifdata_ms_len = acq_samplerate/1000;
conv_acq_data = acq_ifdata;
down_fre = system_fre - ifdata_center_fre;
acq_ifdata_len = acq_ms/1000*acq_samplerate;
end
%% acq process
acq_check_cnt = 1;
for prn = acq_prn_list
ca_code = gnss_code(:,prn);
codeValueIndex = ceil(code_rate*(1:code_len*acq_samplerate/code_rate)/acq_samplerate);
ca_code_resample = ca_code(codeValueIndex); % resample C/A code
ca_code_complex = complex(ca_code_resample,0);
dop_cnt = 1;
isacq_cnt = 0;
for dop = acq_fre_min : acq_fre_step: acq_fre_max
%% set search parameter
temp_noncoh = zeros(acq_ifdata_ms_len*2-1,1);
temp_coh = zeros(acq_ifdata_ms_len*2-1,1);
%%NCO
fncos = dop + down_fre;
dt = 1/acq_samplerate;
t = (0:acq_ifdata_len-1)*dt;
sin_x = sin(2*pi*fncos*t);
cos_x = cos(2*pi*fncos*t);
x_drc = complex(cos_x', sin_x') .* conv_acq_data;
for i = 1 : noncoh_num
for j = 1 : coh_num
ms_cnt = (i - 1)*coh_num + j;
temp_index = 1 + (ms_cnt-1)*acq_ifdata_ms_len : acq_ifdata_ms_len*ms_cnt;
xdr_ms_ij = x_drc(temp_index);
%Normalized results, parallel code correlator
[fft_code_result, x] = xcorr(xdr_ms_ij, ca_code_complex, 'coeff');
%coherent
temp_coh = temp_coh + fft_code_result;
end
%%捕获
temp_noncoh = temp_noncoh + abs(temp_coh);
temp_coh = zeros(acq_ifdata_ms_len*2-1, 1);
end
%% acq check
temp_noncoh_result(:,dop_cnt) = temp_noncoh;
acq_check_result = acq_check(dop, temp_noncoh, acq_snr_thd, code_phase_step);
acq_check_result.dop_cnt = dop_cnt;
if acq_check_result.isacq == 1
disp(['system_type:', num2str(system_type), ' acq prn:', num2str(prn), ' acq_snr:', num2str(acq_check_result.acq_snr), ' codephase: ', num2str(acq_check_result.codephase), ' doppler: ', num2str(acq_check_result.dop)]);
%% save acq_check_list
acq_check_result.system_type = system_type;
acq_check_list{acq_check_cnt, 1}.acq_check_result = acq_check_result;
acq_check_list{acq_check_cnt, 1}.x = x;
acq_check_list{acq_check_cnt, 1}.prn = prn;
%%
isacq_cnt = isacq_cnt + 1;
acq_check_cnt = acq_check_cnt + 1;
acq_result.acq_status = 1;
end
dop_cnt = dop_cnt + 1;
end
%% save noncoh result
if isacq_cnt > 0
for check_cnt = acq_check_cnt - isacq_cnt : acq_check_cnt - 1
acq_check_list{check_cnt, 1}.temp_noncoh_result = temp_noncoh_result;
end
end
end
if acq_result.acq_status == 1
acq_result.acq_check_list = acq_check_list;
end
end
end
%% acq check
function acq_check_result = acq_check(dop, temp_noncoh, acq_snr_thd, code_phase_step)
[~,index_max] = max(temp_noncoh, [], 1); %%按列
index_len = floor(1/code_phase_step/2);
sum1 = sum(temp_noncoh);
sum2 = sum(temp_noncoh(index_max - index_len : index_max + index_len)); %% code step = 0.5
len = length(temp_noncoh(:, 1));
acq_snr = temp_noncoh(index_max)*(len - index_len*2-1)/(sum1 - sum2);
%% acq check
if acq_snr > acq_snr_thd
acq_check_result.isacq = 1;
else
acq_check_result.isacq = 0;
end
acq_check_result.index_max = index_max;
acq_check_result.codephase = code_phase_step*(index_max - floor(index_max/2));
acq_check_result.acq_snr = acq_snr;
acq_check_result.dop = dop;
end
%% resample data
function conv_data = downconvertsample(data, downfre_config)
%% f_result = fc - lco_f
lco_f = downfre_config.fc - downfre_config.fc_down;
ratio_f = downfre_config.fs/downfre_config.fs_down;
%% genertate local signal for convert
len = length(data);
t = 0:1/downfre_config.fs: len*1/downfre_config.fs - 1/downfre_config.fs;
t = t';
real = cos(2*pi*t*lco_f)*32;
imag = sin(2*pi*t*lco_f)*32;
loc = complex(real, imag);
ifdata = complex(data(:,1), data(:, 2));
%% fre down convert
temp_data = loc.*ifdata;
%% down sample rate
temp_ifdata = resample(temp_data, 4, ratio_f*4); %% for 8.184E6
conv_data = temp_ifdata;
end
3 测试结果
3.1 GPS L1CA测试结果
GPS L1CA参数设置：
- 数据采样率2.046 MHz
- 1 ms相干积分，1次相干积分，20次非相干积分
- 捕获步进500Hz
测试结果 2-D图：
测试结果3-D图
3.2 BDS B1I测试结果
BDS B1I参数设置：
- 数据采样率4.092MHz
- 1 ms相干积分，1次相干积分，20次非相干积分
- 捕获步进500Hz
注意：由于BDS 非GEO卫星调制了NH码，无法直接进行大于1ms的相干积分。
测试结果 2-D图：
测试结果3-D图：
4 数据及代码链接
代码和测试数据: gitee link
1042
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
