基于Matlab的TEQC绘图程序代码.docx
《基于Matlab的TEQC绘图程序代码.docx》由会员分享，可在线阅读，更多相关《基于Matlab的TEQC绘图程序代码.docx（24页珍藏版）》请在冰豆网上搜索。
基于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〈length（A）；
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列分别为年月日时分秒。
构造cal时可以省略
升级会员 