基于Matlab的TEQC绘图程序代码.docx

基于Matlab的TEQC绘图程序代码.docx

《基于Matlab的TEQC绘图程序代码.docx》由会员分享，可在线阅读，更多相关《基于Matlab的TEQC绘图程序代码.docx（23页珍藏版）》请在冰豆网上搜索。

基于Matlab的TEQC绘图程序代码

%在matlab下新建一个m文件，将以下代码直接拷贝进去，即可执行。

%需要一个TEQC生成的plot文件作为参数

functionout=teqcplot3（files）;

%读取TEQC生成的Plot文件，绘制数据图表，支持Copmact、Compact2、Compact3格式

%选取一个TEQC的Plot文件

%格式说明

%*.sn1载波L1的信噪比Signaltonoiseratio（S/N）

%*.sn2载波L2的信噪比Signaltonoiseratio（S/N）CarrierL2

%*.iod*.d12*.d21电离层延迟观测值变化率（米/秒）Derivativeofionosphericdelayobservable（m/s）

%*.ion*.i12*.i21电离层延迟观测值（米）Ionosphericdelayobservable（m）

%*.mp1*.m12载波L1的多路径误差MultipathCarrierL1

%*.mp2*.m21载波L2的多路径误差MultipathCarrierL2

%*.azi卫星方位角（°）Satelliteazimuthaldata（degrees）

%*.ele卫星高度角（°）Satelliteelevationdata（degrees）

ifnargin==0

[filen,path]=uigetfile（'*.sn1;*.sn2;*.iod;*.ion;*.mp1;*.mp2;*.azi;*.ele;*.i12;*i21;*.m12;m21;*.d12;*d21',...

'请选择TEQC报告文件：

'）;

else

[path,filen,ext]=fileparts（files）;

path=[path'\'];

filen={[filenext]};

end

%读取这个文件

file=char（filen）;

%按行读取文件至数组A

[A]=importdata（[pathfile],'\t'）;

%定义SAT，存放卫星数据

%GPS有32颗卫星，存放序号1-32，

%GLONASS有32颗卫星，存放序号33-64，

%BEIDOU有35颗卫星，存放序号65-99

SAT（1:

length（A）,1:

99）=NaN;

%sats（1:

length（A）,1:

99）=NaN;

%存放采样时间，单位秒

tsec（1:

length（A））=NaN;

%读取文件的第一行

filelx=A{1};

%判断是哪种格式

switchfilelx

case'COMPACT'

%读取数据采样间隔

t_samp=char（A（3））;

%读取开始时间

mjl=char（A（4））;

%读取数据采样间隔

T_SAMP=str2num（t_samp（max（find（t_samp==''））:

end））;

%读取数据采样开始时间

MJL_START=str2num（mjl（max（find（mjl==''））:

end））;

%转成时间序列数字，dateserialnumber，从0000年1月1日0时0分0秒开始计算的十进制天数

MJD_START=MJL_START+678941.999999741;

%i为行号

n=1;i=5;

case'COMPACT2'

%读取数据采用间隔

t_samp=char（A

（2））;

%读取开始时间

mjl=char（A（3））;

%读取数据采样间隔

T_SAMP=str2num（t_samp（max（find（t_samp==''））:

end））;

%读取数据采样开始时间

MJL_START=str2num（mjl（max（find（mjl==''））:

end））;

%转成时间序列数字，dateserialnumber，从0000年1月1日0时0分0秒开始计算的十进制天数

MJD_START=MJL_START+678941.999999741;

n=1;i=4;

case'COMPACT3'

%读取开始时间

t_start=char（A

（2））;

t_start=deblank（t_start）;

s=splitstr（t_start,'**',6）;

%t_start_time=[char（s{2}）'年'char（s{3}）'月'char（s{4}）'日'char（s{5}）'时'char（s{6}）'分'num2str（str2num（char（s{7}））,'%02d'）'秒'];

%获取采样的开始时间，2013,12，7,03,05,55

t_s_time=[str2num（char（s{2}））,str2num（char（s{3}））,str2num（char（s{4}））,str2num（char（s{5}））,str2num（char（s{6}））,str2num（char（s{7}））];

n=1;i=3;

otherwise

disp（'数据格式存在问题'）;

return

end

%n=1;i=3;

%sats=str2num（A{i}）;

snyggfilen=strrep（filen,'_','-'）;

%生成进度条

h=waitbar（0,['正在读取数据，请稍等……'char（snyggfilen）]）;

whilei

waitbar（i/length（A）,h）;drawnow

%读取卫星编号数据行采样时间卫星数量卫星编号COMAPCT3

satbhstr=char（A（i））;

s=splitstr（satbhstr,'**',32）;

%如果卫星编号数据为0，说明数据缺失，退出本函数

ifstrtrim（satbhstr）=='0'

disp（'数据为空！

'）;

close（h）;

return

end

switchfilelx

case'COMPACT3'

%记录采样时间，单位秒

tsec（n）=str2num（s{1}）;

ifs{2}=='-1'

%如果为-1，使用上次的卫星编号行

satbhstr=oldsatbh;

s=splitstr（satbhstr,'**',32）;

else

oldsatbh=satbhstr;

end

%获取当前卫星数量

satcount=str2num（s{2}）;

case{'COMPACT','COMPACT2'}

%记录采样时间，单位秒

tsec（n）=（n-1）*T_SAMP;

ifs{1}=='-1'

%如果为-1，使用上次的卫星编号行

satbhstr=oldsatbh;

s=splitstr（satbhstr,'**',32）;

else

oldsatbh=satbhstr;

end

%获取当前卫星数量

satcount=str2num（s{1}）;

end

%读取卫星数据行

satdatastr=char（A（i+1））;

sdata=splitstr（satdatastr,'**',32）;

fork=3:

satcount+2

switchfilelx

case'COMPACT'

%如果是COMPACT格式，卫星编号前面没有字母，默认是GPS卫星，在前面增加字符G

satbh=['G'char（s{k-1}）];

case'COMPACT2'

%如果是COMPACT2格式，卫星编号行第2个数据开始是卫星编号

satbh=char（s{k-1}）;

case'COMPACT3'

%如果是COMPACT3格式，卫星编号行第3个数据开始是卫星编号

satbh=char（s{k}）;

end

switchsatbh

（1）;

case'G'

%如果是GPS卫星，获得与编号对应的存储序号，范围1-32

index=str2num（satbh（2:

3））;

case'R'

%如果是GLONASS卫星，获得与编号对应的存储序号，范围33-64

index=32+str2num（satbh（2:

3））;

case'C'

%如果是BEIDOU卫星，获得与编号对应的存储序号，范围65-99

index=64+str2num（satbh（2:

3））;

otherwise

%其他情况，获得与编号对应的存储序号，范围1-32

index=str2num（satbh（2:

3））;

end%switch

%获得对应的卫星数据

sdatastr=char（sdata{k-2}）;

%如果卫星数据不是数值型，用0替代

ifisempty（str2num（sdatastr））==1

sdatastr='0.000';

end

%将卫星数据存入SAT数组的对应位置

SAT（n,index）=str2num（sdatastr）;

end%for

n=n+1;

%根据数据文件结构，一行是卫星编号，下一行就是对应的卫星数据

%程序一次读取2行数据

i=i+2;

end%while

%数据读取完毕，关闭进度条

close（h）;

%获取已使用的卫星编号，返回给sat_bh

sat_bh=getsatbh（SAT）;

%增加最后一个结束字符

sat_bh_temp=[sat_bh'{'End'}];

sat_bh_end=sat_bh_temp';

%删去所有都为NaN值的列，返回给sat_data

sat_data=delnancol（SAT）;

%删去所有NaN值的采样时间列

t_seconds=delnancol（tsec）;

%删去所有都为NaN值的行，返回sat_datas

sat_datas=delnanrow（sat_data）;

%增加一个nan列

[row,col]=size（sat_datas）;

bb（1:

row,1）=NaN;

sat_datas=[sat_datasbb];

%将含有NaN值的数据替换为0，返回给s_data

%s_data=repnan20（sat_data）;

%将采样时间转成顺序日期serialdatenumber

switchfilelx

case{'COMPACT','COMPACT2'}

t_s_jd_time=MJD_START;

case'COMPACT3'

t_s_jd_time=datenum（t_s_time）;

end

fori=1:

length（t_seconds）

%将每次采样时间转换成相应的顺序日期

t_jd_time（i）=t_s_jd_time+t_seconds（i）/60/60/24;

end

%绘制图表++++++++++++++++++++++++++

%根据不同的文件类型，设置坐标轴的范围

[type,maxy,miny]=get_filetype（file）;

figure;boxon;holdon

[row,col]=size（sat_datas）;

%绘制渐变彩色图

pcolor（t_jd_time',[1:

col],sat_datas'）;

set（gcf,'renderer','zbuffer'）;

shadingflat

set（gca,'xticklabel',[t_jd_time]）;

set（gca,'xlim',[t_jd_time

（1）t_jd_time（end）]）

%t_start_times=[t_start_h':

't_start_m':

't_start_s];

dateaxis（'X',13）

cbar（'v',[minymaxy],type）;

set（gca,'ylim',[1col+1]）

set（gca,'ytick',[1.5:

1:

col+0.5]）

set（gca,'yticklabel',sat_bh_end）;

set（gca,'fontsize',7）;

colormap（flipud（jet））;

caxis（[minymaxy]）;

%t_end_time=cal2et（t_s_time,t_seconds（end））;

xlabel（[datestr（t_jd_time

（1））'|--------采样间隔:

'num2str（t_seconds

（2）-t_seconds

（1））'秒--------|'datestr（t_jd_time（end））]）

ylabel（'卫星编号'）

%timestr=secs2hms（length（sat））;

T=title（['TEQC报告文件:

'strrep（file,'_','-'）]）;

set（T,'fontsize',8）

%out.（file（end-2:

end））=sat;

%out.T_samp=T_SAMP;

out.Start=datestr（t_jd_time

（1））;

out.Stop=datestr（t_jd_time（end））;

%+++++++++++++++++++++++++++++++++++++++

%以下是本函数中所使用的一些相关函数

%++++++++++++++++++++++++++++++++++++++++++++++++++

%+++++++++++++++++++++++++++++++++++++++++++++++++

functions=splitstr（str,split,num）;

%用指定分隔符分隔字符串函数

%str='name02010-05-0633QtBKB';

%'**'表示特殊分隔符回车、空格、换行、制表符等

%num表示分隔字符串的个数

string=cell（num,1）;

index=1;

temp=[];

ifsplit=='**'

%特殊分隔符

fori=1:

length（str）

ifisspace（str（i））

ifisempty（temp）==false

string{index}=temp;

index=index+1;

temp=[];

end

else

temp=[temp,str（i）];

end

end

else

%常规分隔符

fori=1:

length（str）

ifstr（i）==split

ifisempty（temp）==false

string{index}=temp;

index=index+1;

temp=[];

end

else

temp=[temp,str（i）];

end

end

end

ifisempty（temp）==false

string{index}=temp;

end

s=string;

%+++++++++++++++++++++++++++++++++++

functionB=delnancol（A）;

%删除矩阵中所有值都是nan的列

[xy]=size（A）;

fori=y:

-1:

1

ifall（isnan（A（:

i）））==1

A（:

i）=[];

end

end

B=A;

%+++++++++++++++++++++++++++++++++++

functionB=delnanrow（A）;

%删除矩阵中所有值都是nan的行

[xy]=size（A）;

fori=x:

-1:

1

ifall（isnan（A（i,:

）））==1

A（i,:

）=[];

end

end

B=A;

%+++++++++++++++++++++++++++++++++++

functionB=repnan20（data）;

%将矩阵中的NAN值用0替代

[datas,features]=size（data）;

fork=1:

features

fori=1:

datas

ifisnan（data（i,k））==1

data（i,k）=0;

end

end

end

B=data;

%++++++++++++++++++++++++++++++++

functionsatbh=getsatbh（satdata）;

%获取卫星数据中的卫星编号，返回一个字符串数组

s=cell（32,1）;

temp=[];

index=1;

fori=1:

99

ifall（isnan（satdata（:

i）））~=1

switchi>32

case0

temp=['G',num2str（i,'%02d'）];

case1

ifi>64

temp=['C',num2str（i-64,'%02d'）];

else

temp=['R',num2str（i-32,'%02d'）];

end

end

s{index}=temp;

index=index+1;

end

end

ifindex>1

satbh=s（1:

index-1）;

else

satbh=[];

end

%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%FUNCTION:

GETFILETYPEINFO+++++++++++++++++++++++++++++++++++++

function[out,maxy,miny]=get_filetype（teqfile）;

%根据文件类型，设置坐标轴的范围

[path,name,ext,ver]=fileparts（teqfile）;

switchext

case'.sn1'

%out='SignaltonoiseratioS/NL1';

out='L1载波上的信噪比（dBHz）';

maxy=60;

miny=15;

case'.sn2'

%out='SignaltonoiseratioS/NL2';

out='L2载波上的信噪比（dBHz）';

maxy=60;

miny=15;

case{'.mp1','.m12'}

%out='MultipathL1';

out='L1多路径观测值';

maxy=1;

miny=-1;

case{'.mp2','.m21'}

%out='MultipathL2';

out='L2多路径观测值';

maxy=1;

miny=-1;

case{'.iod','.d12','d21'}

%out='Derivativeofionosphericdelayobservable（m/s）';

out='电离层延迟观测值变化率（米/秒）';

maxy=1;

miny=-1;

case{'.ion','.i12','i21'}

maxy=2;

miny=-2;

%out='Ionosphericdelayobservable（m）';

out='电离层延迟观测值（米）';

case'.ele'

maxy=90;

miny=0;

%out='Satelliteelevationdata';

out='卫星高度角（°）';

case'.azi'

maxy=180;

miny=-180;

%out='Satelliteazimuthaldata';

out='卫星方位角（°）';

otherwise

disp（'Somethingswrong..!

'）

end

%FUNCTION:

PLACEAMODIFIEDCOLORBAR++++++++++++++++++++++++++++++

functionCB=cbar（loc,range,label）;

%.............................................................

%CB=cbar（loc,range,label）

%placesacolorbarat:

%loc='v'inverticalor'h'inhorizontal

%positionincurrentfigurescaledbetween:

%range=[minmax]witha:

%label='string'.

%

%fontsizeisreducedto10andwidthofbarishalfdefault.

%

%Example:

[X,Y,Z]=peaks（25）;

%range=[min（min（Z））max（max（Z））];

%pcolor（X,Y,Z）;

%cbar（'v',range,'Elevation（m）'）

%.............................................................

caxis（[range

（1）range

（2）]）;

switchloc

case'v'

CB=colorbar（'vertical'）;

set（CB,'ylim',[range

（1）range

（2）]）;

POS=get（CB,'position'）;

set（CB,'position',[POS

（1）POS

（2）0.03POS（4）]）;

set（CB,'fontsize',8）;

set（get（CB,'ylabel'）,'string',label）;

set（CB,'box','on'）

case'h'

CB=colorbar（'horizontal'）;

set（CB,'xlim',[range

（1）range

（2）]）;

POS=get（CB,'position'）;

set（CB,'position',[POS

（1）POS

（2）POS（3）0.03]）;

set（CB,'fontsize',8）;

set（get（CB,'xlabel'）,'string',label）

end

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

functionjd=cal2jd（cal）

%cal2jd1将公历年月日时分秒转换到儒略日。

%jd=cal2jd（cal）返回儒略日

%cal：

1x6矩阵，6列分