clear all;clc;cla;close all;
format longg; %%15位有效数
%% 一、 设定GPS数据文件参数 %%%%%%%%%%%%%%
gps_bd = 0; % 0:gps, 1:bd
n_mean_count = 200; %bias count
if gps_bd == 0
% 这个中频数据是 bit8
filename = '.\data\GPS_sim_5mfs_if1.42m.txt';% 相对路径 当前工程文件夹下的data文件夹下的GPS_sim_5mfs_if1.42m.txt文件
fs = 5e6; % 采样率★书本260页实际采样在5E6
fc = 1.42e6;% 采样数据中频★
sv = 1:32;% 要捕获跟踪的卫星号 若不知道有哪些星可以设置为 sv=1:32 但这样处理速度会比较慢★
Readtype = 'bit8';%数据读取位宽，要注意 不同的中频信号采样位是否一致，有符号采样时，高位为符号位，低位为幅值位
% filename = 'D:\9.WorkInTheLab\0.Project\0.GNSS_SignalQualityMonitor\1.data\2771\data\GPSL1_8bit_fs16.369M_fc3.996M.bin';
% fc = 3.996e6;
% fs = 16.369e6;
% sv = 28;
% Readtype = 'bit8';
% filename = 'C:\MATLAB7\work\BDSB1-8bit-PwrSpectAbnorm-202006172224.bin';
% fc = 3.996e6;
% fs = 16.369e6;
% sv = 1:32;
% Readtype = 'bit8';
disp 'gps acq'
else
filename = 'C:\MATLAB7\work\usbdata.bin';
fc = 3.999875e6;
fs = 16.369e6;
sv = 1:32;
Readtype = 'bit8';
disp 'bd acq'
end
%% 二、 任务配置 %%%%%%%%%%%%%%
task_AcqOnly= 0 ; %捕获标志，为1时只进行捕获，为0时则进行捕获跟踪
task_Acq_Track=1; %捕获跟踪标志，为1 － 则进行捕获、跟踪；为 0 则调用以前结果进行NAV定位
track_acqresultsName='capture_test.mat';%%捕获跟踪后生成的数据保存名字,task_Acq_Track = 0 时则直接调用此数据进行定位解算。注意每次用不同的中频数据时取不同的名字，以免覆盖。
task_frameFindChk= 1 ; %找子帧标志， 1则找子帧，其他则不进行 0
task_nav = 1 ; %定位解算标志 1则参与定位解算，其他则不进行 0
plotAcq_AND_Channel= 1 ; %捕获通道画图标志，1则画图，其他则不进行
plot_Acq_ACF= 1 ; %捕获三维图和相关累加画图标志，1则画图，其他则不进行
if task_Acq_Track==1 %捕获跟踪标志，为1 － 则进行捕获、跟踪；为 0 则调用以前结果进行NAV
task_load_results= 0 ;
else
task_load_results= 0 ;
end
%% 三、 捕获、跟踪参数设置 %%%%%%%%%%%%%%
% file_sec=0.400; %跟踪时间★，单位为秒，理论上大于30秒后方可能定位解算，一般取33来验证，也可长时间跑看定位效果 %%!!!!!!注意如果有时数据长度不够的话，会报错。如此时的数据GPS_sim_5mfs_if1.42m.txt用36s就会报错，用33则不会。
file_sec=33; %跟踪时间要这么长？★
channels=8; %通道数设置★8-12书上，但后面真正用到的不一定是8,可能<8
acq_th=2.5; %捕获门限★
K=1; %非相干累加次数★
fc_lo=fc-5000;fc_hi=fc+5000;step=100;%doppler bin★
nav_interval=0.001; %导航定位解算时间间隔，单位为秒★
ts=1/fs; %采样时间★
n=fs/1000; n_mli=n; %1ms的采样点数 ★
nn=0:n-1; %采样点范围 ★
num_sec=0.011; %读取文件中的长度设置，★
n_samples=fs*num_sec; %读取文件中的采样点个数，如当前设置为10ms采样数据
plot_acq_bar_par = []; %捕获柱状画图卫星SNR计数
%% 四、 读取 真实采集数据或调用以前的处理结果%%%%%%%%%%%%%%%%%%%%%
if task_Acq_Track==1
fid = fopen(filename,'rb'); %打开filename路径读取文件
skipNumberOfBytes = fs * 3; %为什么是3？
fseek(fid,skipNumberOfBytes, 'bof');
data_sum=fread(fid,n_samples,Readtype); %以Readtype为文件格式，以n_samples为 数据长度已读取文件保存到data中
fclose(fid); %关闭文件
data_sum=data_sum';%%%对读取数据进行处理，目的是把数据的值缩小一半缩小计算量，也可不处理拿来直接用
n=4;
data=data_sum;
[data_bias,mean_data] = bias_control(data,n_mean_count);
data = data_bias;
plotPSD(data,fs);grid on; %功率画图
xlabel('频率/MHz'); ylabel('功率谱密度/dB·Hz');
title ('信号功率谱密度');
end
%%%%%%%%*********.load 以前的结果进行定位解算. .**********%%%%%%%%
if task_load_results==1
load(track_acqresultsName); %导入捕获跟踪存储的数据直接用来导航定位解算
end
%% 五、 ★捕获%%%%%%%%%%%%%%%%%%%%
if task_Acq_Track==1
disp '1..Acqisiting.!----'
acq=[]; %acq初始化
j=1; %捕获卫星个数变量，用于保存捕获通道参数指向路径
fprintf('(');
for i=1:length(sv) %捕获设定的卫星编号
svnum=sv(i); %当前捕获的卫星编号
if gps_bd == 0
% 输入参数：中频数据，卫星编号，采样频率，频率搜索起点，频率搜索步长，频率搜索截止点，非相干累加次数，捕获门限
% 输出参数：捕获成功标志，捕获的载波频率，捕获的码相位，非相干累加值，非相干累加值峰均比
% [acqOK,f,ph,acf,PAR]=GPSacquisition(data,svnum,fs,fc_lo,step,fc_hi,K,acq_th);
[acqOK,f,ph,acf,PAR] = GPSacquisition_s(data,svnum,fs,fc_lo,step,fc_hi,K,acq_th);
else
[acqOK,f,ph,acf,PAR]=BDacquisition(data,svnum,fs,fc_lo,step,fc_hi,K,acq_th);
end
plot_acq_bar_par(i)=PAR;
% ★捕获成功时，缓存捕获信息：非相干累加值，卫星编号，牵引载波频率，码相位
% 画捕获三维图和非相干累加值图
if acqOK
figure;
acfhandler=mesh(acf);set(acfhandler,'EdgeColor','none','Facecolor','interp');
text=sprintf('SV %i',svnum);title(text); xlabel('码相位');ylabel('载波多普勒');
[m fbin]=max(max(acf'));% 行号即是载波多普勒
[m ccn]=max(max(acf));% 列号即是码相位
% figure; plot(abs(acf(:,ccn)),'-o');
% xlabel('多普勒');ylabel('非相干积分值');
% figure; plot(abs(acf(fbin,:)),'-o');
% xlabel('码相位');ylabel('非相干积分值');
% [f_pull,dfreq]=gps_pvt_fpullin(data,svnum,f,fs,ph); % 5hz 高分辨率
acq(j,1)=PAR;
acq(j,2)=sv(i); %卫星编号
acq(j,3)=f;
acq(j,4)=ph; %码相位
j=j+1; %每捕获一颗星，通道号+1
fprintf('%d，', sv(i)); %捕获成功是则输出卫星号
else
fprintf('.'); %未捕获成功时则输出"."
end
end
fprintf(')\n');
if isempty(acq)
disp '---acq is failed---';
return;
else
plot_acq_bar(plot_acq_bar_par,acq_th); %捕获柱状画图
disp '2..Configure channels for track: '
disp ' _______________________________________________________'
disp '| CH | PRN | Frequency | Doppler | CodePhase | PAR |';
channel_num=min(channels,size(acq,1)); %取决于初始通道数和捕获卫星数
if plotAcq_AND_Channel==1
for i=1 : channel_num
fprintf('| %2d | %2d | %6.2f | %6.0f | % 6d | %2.2f | \n',...
i, ... %显示通道号
acq((i),2), ... %显示卫星编号
acq((i),3), ... %载波频率
acq((i),3) - fc, ... %多普勒频移
acq((i),4),... %码相位
acq((i),1)); %非相干累加幅值
end
end
disp ' --end of Acqisition!AND Sort acqresults by PAR-----';
[acq_SNRdescend]=gps_SoftReceiver_sort_acqdescending(acq);%把捕获的结果按SNR的顺序由大到小重新排列
acq=acq_SNRdescend;
end
end
%% 六、★跟踪并处理、保存结果
if task_AcqOnly==0 %到本脚本的结束位置
if (task_Acq_Track==1&&task_AcqOnly==0)
%%6.1---跟踪初始化 % %
disp '3..Tracking.!----'
length_ms=file_sec*1000; %跟踪时间长度
Ips=zeros(channel_num,length_ms); %I路信号初始化
Qps=zeros(channel_num,length_ms); %Q路信号初始化
code_shifts=zeros(channel_num,length_ms); %码相位偏移初始化
navidatas=zeros(channel_num,length_ms/20+2); %电文数据比特长度初始化
navdata_indexs=zeros(channel_num,length_ms/20+2); %电文数据比特长度深度初始化
carr_Freqs=zeros(channel_num,length_ms); %载波频率初始化
carr_remPhases=zeros(channel_num,length_ms); %载波角频率初始化
code_Freqs=zeros(channel_num,length_ms); %码频率初始化
code_remPhases=zeros(channel_num,length_ms); %码角频率初始化
%%6.2******跟踪－输出导航电文 + 数据解调＋1－1 并主意这里的跟踪相位*%%%%%%%%%
Bit_sync_OK =0;
Bit_sync =zeros(channel_num,length_ms);
for i=1:channel_num
fprintf(' %d ',i) %输出当前跟踪通道
phase=acq((i),4); %码相位
%函数：gps_pvt_track_beijingNCO 跟踪环
% 输入：[路径文件名，卫星号，采样频率，捕获载波频率，捕获码相位，跟踪时间，读取数据位宽，捕获通道总数，当前通道号]
% 输出：[I路数据，Q路数据，码相位偏移地址，EPL功率值，包含多普勒频移载波频率，载波角相位，码频率，码角相位]
[Ip,Qp,code_ph_shift,acf_EPL,carr_Freq,carr_remPhase,code_Freq,code_remPhase]=gpsTracking(filename,acq((i),2),fs,acq((i),3),phase,file_sec,Readtype,skipNumberOfBytes,channel_num,i);
% Ips(i,1:length(Ip))=Ip; %Ip路数据
% Qps(i,1:length(Qp))=Qp; %Qp路数据
% code_shifts(i,1:length(code_ph_shift))=code_ph_shift; %码相位偏移地址
% carr_Freqs(i,1:length(carr_Freq))=carr_Freq; %包含多普勒频移的载波频率
% carr_remPhases(i,1:length(carr_remPhase))=carr_remPhase; %载波角相位
% code_Freqs(i,1:length(code_Freq))=code_Freq; %码NCO
% code_remPhases(i,1:length(code_remPhase))=code_remPhase; %码角相位
end
% save (track_acqresultsName, 'channel_num', 'navidatas','code_shifts', 'Qps', 'Ips', ...
% 'navdata_indexs','Bit_sync','acf_EPL' ,'acq','fs','carr_Freqs','carr_remPhases','code_Freqs','code_remPhases','fc');%注意每次保存后更改data的名字 以免覆盖
% fprintf('\n');
disp '..End of Track AND track_acqresults saved!----'
else
disp 'task_track=0'
end
figure; subplot(2,1,1),plot(Ips,'-o');
subplot(2,1,2), plot(Qps,'-*');
end
%
版权声明：本文为LWN7321原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。