讀取GPS觀測數(shù)據(jù)O文件地matlab編程

functionHeadO=ReadObsHead[fname,fpath]=uigetfile('*.*O','選擇一個O文件');HeadO.O_filename=strcat(fpath,fname);fid=fopen(HeadO.O_filename,'rt');if(fid==-1)msgbox('文件不存在','警告','warn');return;endk=0;while(1)s=fgets(fid);if(strfind(s,'RINEXVERSION'))HeadO.Rtype=str2double(s(6:10));elseif(strfind(s,'APPROXPOSITIONXYZ'))HeadO.XYZ(1)=str2double(s(1:14));HeadO.XYZ(2)=str2double(s(15:28));HeadO.XYZ(3)=str2double(s(29:42));elseif(strfind(s,'MARKERNAME'))HeadO.staname=s(1:4);elseif(strfind(s,'MARKERNUMBER'))HeadO.stanum=s(1:4);elseif(strfind(s,'#/TYPESOFOBSERV'))HeadO.PRNtype='G';HeadO.nObstype=str2double(s(6));fori=1:HeadO.nObstypeHeadO.Obs{i}={s(i*6+5:i*6+6)};endelseif(strfind(s,'SYS/#/OBSTYPES'))k=k+1;HeadO.PRNtype(k)=char(s(1));HeadO.nObstype(k)=str2double(s(5:6));if(HeadO.nObstype(k)>13)fori=1:13HeadO.Obs{k,i}=s(i*4+4:i*4+6);ends=fgets(fid);fori=1:HeadO.nObstype(k)-13HeadO.Obs{k,i+13}=s(i*4+4:i*4+6);endelsefori=1:HeadO.nObstype(k)HeadO.Obs{k,i}=s(i*4+4:i*4+6);endendelseif(strfind(s,'ENDOFHEADER'))break;endend下面為讀取數(shù)據(jù)程序，供參考

%ReadRINEXObservationData

nepoch=0;

nobs=0;

nlli=0;

rnx.epoch=[];

rnx.obs=[];

rnx.lli=[];

while1

iffeof(fid)==1,break,end

line=fgetl(fid);

ifisempty(deblank(line)),break,end

ifversion<3

%Version2RINEX

%getepochflagandnumberofsatellites

EpStr=line(1:26);

EpFlg=sscanf(line(27:29),'%d');

NumSat=sscanf(line(30:32),'%d');

ifEpFlg>1

%epochflags

%

0-o.k.

%

1-powerfailurebetweenpreviousandcurrentepoch

%eventflags(epochflagcontinued)

%

2-startmovingantenna

%

3-newsiteoccupation(endofkinematicdata);atleastarecord

%

withanewmarkernamemustfollow

%

4-newheaderinfofollows

%

5-externalevent(epochissignificant)

%

6-cyclesliprecordsfollow(sameformatasobservations)

%

%thenumberofrecordswhichfollowisgivenbyNumSat

fori=1:NumSat

line=fgetl(fid);

end

continue

end

%getepochtime

t=sscanf(EpStr,'%f',6)';

ift(1)<80

t(1)=2000+t(1);

else

t(1)=1900+t(1);

end

%keyboard

%-adaptionrvbree25112008---

%

date

=datenum(t(1),t(2),t(3));%Theserialdatenumbersfrom1-Jan-0000tothedayofobservation

%

gpsweek

=fix((date-datenum(1980,1,6))/7);

%

secsinweek

=(date-datenum(1980,1,6)-gpsweek*7)*24*3600+t(4)*3600+t(5)*60+t(6);

%

%

EpTime

=secsinweek;

%

GPSW

=gpsweek;

%EpTime=datenum(t(1),t(2),t(3),t(4),t(5),t(6));

%----endadaption-----------------------------------------

%EpTime=datenum(t(1),t(2),t(3),t(4),t(5),t(6));

%readoptionalclockerror***wedonotneedtheclockerror***

%ClkOffset=sscanf(line(69:end),'%f');

%getsatelliteid's

%-sprncontainsthe3lettercodesforthesatellites

%-prncontainstheuniversalsatelliteid's

sprn=line(33:min(end,68));

forj=2:ceil(NumSat/12)

line2=fgetl(fid);

sprn=[sprnline2(33:min(end,68))];

end

prn=rxsatnum(sprn);

%gettheobservations,lliandsnrforallsatellites

obs=zeros(NumSat,NumRxObs);

fori=1:NumSat

line=fgetl(fid);line=[lineblanks(80-length(line))];

forj=2:ceil(NumRxObs/5)

line2=fgetl(fid);

line=[lineline2blanks(80-length(line2))];

end

%gettheobservations

p=0;

forj=1:NumRxObs

obsi=sscanf(line(p+1:p+14),'%f');

%testformissingobservations

ifisempty(obsi)|obsi==0,obsi=NaN;,end

obs(i,j)=obsi;

p=p+16;

end

%getthelli(lossoflockindicator)forcarrierphasedata

forjj=1:length(icarrier)

j=icarrier(jj);

flg=double(line((j-1)*16+15))-48;

ifrem(flg,2)