rinex观测文件数据读取

【实例简介】读取rinex观测文件中伪距，载波，多普勒测量值
【实例截图】cebr1050.16o
【核心代码】% function: read the body of rinex obs file
clear all
close all
clc

fid=fopen('cebr1050.16o');
if fid==-1
    disp('fail to open obs file!');
end

% 创建结构体obs，使用Map进行映射，prn为键(char)，sv为值(table)
obs=containers.Map;
record_counter=0;

% 读取body记录
while feof(fid)==0
    line=fgets(fid);
    if strcmp(line(1),'>')
        sv_sum=str2double(line(34:35));
        time=line(3:29);
        if sv_sum>0
            for i=1:sv_sum
                line=fgets(fid);
                % 键 -> prn
                key=line(1:3);
                % 值 -> sv(pr cp dop)
                sv.pr=str2double(line(6:17));
                cp=line(21:33);
                if isnan(cp)
                    sv.cp=0;
                else
                    sv.cp=str2double(line(21:33));
                end
                sv.dop=str2double(line(41:49));
                % 值逻辑判断
                if isKey(obs,key)
                    jump_sec=record_counter-length(obs(key));
                    if jump_sec~=0 % 解决掉星后重新锁定的问题
                        temp.pr=0;
                        temp.cp=0;
                        temp.dop=0;
                        tmp_strcut=repmat(temp,1,jump_sec);
                        value=[obs(key) tmp_strcut sv];
                    else % 卫星一直锁定，直接追加
                        %if strcmp(key,'C18')==1
                         %   if length(obs(key))>5713
                          %      disp(time);
                          %  end
                       % end
                        value=[repmat(obs(key),1) repmat(sv,1)]; % mark
                    end
                else % 卫星首次锁定
                    if record_counter~=0 % 第一秒未锁定
                        temp.pr=0;
                        temp.cp=0;
                        temp.dop=0;
                        tmp_strcut=repmat(temp,1,record_counter);
                        value=[tmp_strcut sv];
                    else % 第一秒锁定
                        value=sv;
                    end
                end
                obs(key)=value;
            end
        end
        record_counter=record_counter 1;
    end
end
disp('create struct obs(body) successfully!');

% 绘图
sv_sum=obs.Count; % 可视卫星数量
keySet=keys(obs); % 可视卫星名称
disp(keySet);

v=299792458.0;
bf1=1.561098E9; % bds L1 频率
gf1=1.57542E9; % gps L1 频率
gps_sum=0;bds_sum=0;

for i=1:sv_sum
    prn=char(keySet(i));
    if strcmp(prn(1),'G')==1 % GPS
        gps_sum=gps_sum 1;
    else if strcmp(prn(1),'C')==1 % BDS
            bds_sum=bds_sum 1;
        end
    end
end
subfig1=ceil(sqrt(bds_sum));
subfig2=ceil(sqrt(gps_sum));

x=1:record_counter;
for i=1:sv_sum
    prn=char(keySet(i));
    info=obs(prn);
    pr=zeros(1,record_counter);cp=pr;
    for j=1:length(info)
        pr(j)=info(j).pr;
        cp(j)=info(j).cp;
    end
    if i==1 || i==bds_sum 1
        figure;
    end
    subplot(4,4,1);
    plot(x,pr,'r',x,cp,'g'); 
    title(char(prn));
end