资源描述
基于Matlab的TEQC绘图程序代码【实用文档】doc
文档可直接使用可编辑,欢迎下载
%在 matlab下新建一个m文件,将以下代码直接拷贝进去,即可执行。
%需要一个TEQC生成的plot文件作为参数
function out=teqcplot3(files);
%读取TEQC生成的Plot文件,绘制数据图表,支持Copmact、Compact2、Compact3格式
%选取一个TEQC的Plot文件
% 格式说明
% *。sn1 载波 L1的信噪比 Signal to noise ratio (S/N)
% *。sn2 载波 L2的信噪比 Signal to noise ratio (S/N) Carrier L2
% *.iod *.d12 *.d21 电离层延迟观测值变化率(米/秒) Derivative of ionospheric delay observable (m/s)
% *。ion *。i12 *.i21 电离层延迟观测值(米)Ionospheric delay observable (m)
% *。mp1 *。m12 载波L1的多路径误差 Multipath Carrier L1
% *。mp2 *.m21 载波L2的多路径误差 Multipath Carrier L2
% *.azi 卫星方位角(°) Satellite azimuthal data (degrees)
% *.ele 卫星高度角(°) Satellite elevation data (degrees)
if nargin==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={[filen ext]};
end
%读取这个文件
file=char(filen);
%按行读取文件至数组A
[A]=importdata([path file],’\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};
%判断是哪种格式
switch filelx
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));
%转成时间序列数字,date serial number,从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));
%转成时间序列数字,date serial number,从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)]);
while i〈length(A);
waitbar(i/length(A),h);drawnow
%读取卫星编号数据行 采样时间 卫星数量 卫星编号 COMAPCT3
satbhstr=char(A(i));
s = splitstr(satbhstr,'**',32);
%如果卫星编号数据为0,说明数据缺失,退出本函数
if strtrim(satbhstr)==’0’
disp ('数据为空!');
close(h);
return
end
switch filelx
case 'COMPACT3'
%记录采样时间,单位秒
tsec(n)=str2num(s{1});
if s{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;
if s{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);
for k=3:satcount+2
switch filelx
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
switch satbh(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替代
if isempty(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_datas bb];
%将含有NaN值的数据替换为0,返回给s_data
%s_data=repnan20(sat_data);
%将采样时间转成顺序日期serial date number
switch filelx
case {'COMPACT’,'COMPACT2’}
t_s_jd_time=MJD_START;
case 'COMPACT3'
t_s_jd_time=datenum(t_s_time);
end
for i=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;box on;hold on
[row,col]=size(sat_datas);
%绘制渐变彩色图
pcolor(t_jd_time',[1:col],sat_datas’);
set(gcf,'renderer','zbuffer');
shading flat
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',[miny maxy],type);
set(gca,'ylim',[1 col+1])
set(gca,’ytick',[1。5:1:col+0.5])
set(gca,'yticklabel',sat_bh_end);
set(gca,'fontsize',7);
colormap(flipud(jet));
caxis([miny maxy]);
%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));
%+++++++++++++++++++++++++++++++++++++++
%以下是本函数中所使用的一些相关函数
%++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++++++++++++++++++++++
function s=splitstr(str,split,num);
%用指定分隔符分隔字符串函数
%str='name 0 2010-05-06 33 Qt BKB';
%'**'表示特殊分隔符回车、空格、换行、制表符等
%num表示分隔字符串的个数
string=cell(num,1);
index=1;
temp=[];
if split==’**'
%特殊分隔符
for i=1:length(str)
if isspace(str(i))
if isempty(temp)==false
string{index}=temp;
index=index+1;
temp=[];
end
else
temp=[temp,str(i)];
end
end
else
%常规分隔符
for i=1:length(str)
if str(i)==split
if isempty(temp)==false
string{index}=temp;
index=index+1;
temp=[];
end
else
temp=[temp,str(i)];
end
end
end
if isempty(temp)==false
string{index}=temp;
end
s=string;
%+++++++++++++++++++++++++++++++++++
function B=delnancol(A);
%删除矩阵中所有值都是nan的列
[x y]=size(A);
for i=y:-1:1
if all(isnan(A(:,i)))==1
A(:,i)=[];
end
end
B=A;
%+++++++++++++++++++++++++++++++++++
function B=delnanrow(A);
%删除矩阵中所有值都是nan的行
[x y]=size(A);
for i=x:-1:1
if all(isnan(A(i,:)))==1
A(i,:)=[];
end
end
B=A;
%+++++++++++++++++++++++++++++++++++
function B=repnan20(data);
%将矩阵中的NAN值用0替代
[datas,features]=size(data);
for k=1:features
for i=1:datas
if isnan(data(i,k))==1
data(i,k)=0;
end
end
end
B=data;
%++++++++++++++++++++++++++++++++
function satbh=getsatbh(satdata);
%获取卫星数据中的卫星编号,返回一个字符串数组
s=cell(32,1);
temp=[];
index=1;
for i=1:99
if all(isnan(satdata(:,i)))~=1
switch i>32
case 0
temp=['G',num2str(i,’%02d')];
case 1
if i>64
temp=[’C',num2str(i-64,’%02d’)];
else
temp=['R’,num2str(i—32,'%02d')];
end
end
s{index}=temp;
index=index+1;
end
end
if index〉1
satbh=s(1:index-1);
else
satbh=[];
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% FUNCTION: GET FILETYPE INFO +++++++++++++++++++++++++++++++++++++
function [out,maxy,miny]=get_filetype(teqfile);
%根据文件类型,设置坐标轴的范围
[path,name,ext,ver]=fileparts(teqfile);
switch ext
case ’。sn1'
%out='Signal to noise ratio S/N L1';
out=’L1载波上的信噪比(dBHz)';
maxy=60;
miny=15;
case '。sn2'
%out='Signal to noise ratio S/N L2’;
out='L2载波上的信噪比(dBHz)’;
maxy=60;
miny=15;
case {’.mp1','。m12'}
%out=’Multipath L1’;
out=’L1多路径观测值';
maxy=1;
miny=-1;
case {’.mp2’,’。m21'}
%out='Multipath L2';
out=’L2多路径观测值';
maxy=1;
miny=—1;
case {'.iod’,'.d12','d21'}
%out='Derivative of ionospheric delay observable (m/s)';
out=’电离层延迟观测值变化率(米/秒)’;
maxy=1;
miny=-1;
case {'。ion’,'.i12’,'i21'}
maxy=2;
miny=-2;
%out=’Ionospheric delay observable (m)';
out=’电离层延迟观测值(米)’;
case '。ele’
maxy=90;
miny=0;
%out='Satellite elevation data’;
out=’卫星高度角(°)';
case ’.azi'
maxy=180;
miny=-180;
%out=’Satellite azimuthal data';
out='卫星方位角(°)';
otherwise
disp('Somethings wrong..!')
end
% FUNCTION: PLACE A MODIFIED COLORBAR ++++++++++++++++++++++++++++++
function CB=cbar(loc,range,label);
% 。。....。.。.。.。.。。....。.。..。...。..。.。.。.。.。...。..。。..。.。。..。...
% CB = cbar(loc,range,label)
% places a colorbar at:
% loc = 'v' in vertical or 'h' in horizontal
% position in current figure scaled between:
% range = [min max] with a:
% label = ’string'。
%
% fontsize is reduced to 10 and width of bar is half default.
%
% 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)]);
switch loc
case 'v'
CB=colorbar('vertical’);
set(CB,’ylim',[range(1) range(2)]);
POS=get(CB,'position’);
set(CB,'position',[POS(1) POS(2) 0。03 POS(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
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function jd=cal2jd(cal)
% cal2jd1ﻩ将公历年月日时分秒转换到儒略日。
% jd=cal2jd(cal) 返回儒略日
% cal:1x6矩阵,6列分别为年月日时分秒。构造cal时可以省略末尾的0
%
% 公元1582年10月4日24:00点之前使用儒略历,公元1582年10月15日00:00点之后使用公历
if length(cal) 〈 6
cal(6)=0;
end
year=cal(1);
month=cal(2);
day=cal(3)+(cal(4)*3600+cal(5)*60+cal(6))/86400;
y = year + 4800; %4801 B.C. is a century year and also a leap year.
if( year 〈 0 )
y =y+ 1; ﻩ% Please note that there is no year 0 A。D。
end
m=month;
if( m 〈= 2 ) % January and February come after December。
m = m+12;
ﻩy = y — 1;
end
e=floor(30.6 * (m+1));
a=floor(y/100); % number of centuries
% 教皇格雷戈里十三世于1582年2月24日以教皇训令颁布,将1582年10月5日至14抹掉。1582年10月4日过完后第二天是10月15日
if( year <1582 )|(year==1582&month<10)|(year==1582& month==10 &day〈15)
b = -38;
else
ﻩb = floor((a/4) - a); % number of century years that are not leap years
end
c=floor(365.25* y); % Julian calendar years and leap years
jd= b + c + e + day - 32167.5;
%+++++++++++++++++++++++++++++++++++++++++++++++++
function cal=jd2cal(J)
% 从儒略日计算公历年月日时分秒
% cal=jd2cal(J)
% 返回的cal是1x6矩阵,6列分别为年月日时分秒
%
% 公元1582年10月4日24:00点之前使用儒略历,公元1582年10月15日00:00点之后使用公历
if (J 〈 1721423.5) % 公元1月1日0时
BC = 1;
else
BC = 0;
end
%start from Julian March 1, 4801 B。C.
if( J < 2299160.5 ) % before 1582.10.4。 24:00 is Julian calender
ﻩj0=floor(J+0.5);
ﻩdd=J+0.5—j0;
else % after 1582。10.15。 00:00 is Gregorian calender
%number of certury years that are not leap year
n1=floor((J—2342031.5)/36524。25/4)+1; %1700。3。1.0
n2=floor((J—2378555.5)/36524.25/4)+1; %1800.3。1。0
n3=floor((J—2415079.5)/36524.25/4)+1; %1900。3。1.0
ﻩj0=n1+n2+n3+J+10;
dd=j0+0。5—floor(j0+0。5);
j0=floor(j0+0.5);
end
j0=j0+32083;
year0=ceil(j0/365.25)—1;
year=year0-4800;
day=j0—floor(year0*365.25);
month=floor((day-0.6)/30.6)+3;
day=day-round((month—3)*30.6);
if month〉12
ﻩmonth=month—12;
year=year+1;
end
year=year—BC;
sec=round(dd*86400);
hour=floor(sec/3600);
sec=sec—hour*3600;
min=floor(sec/60);
sec=sec-min*60;
cal=[year, month, day, hour, min, sec];
%+++++++++++++++++++++++++++++++++++++++++++++++++
function mjd=cal2mjd(cal)
% cal2mjdﻩ将公历年月日时分秒转换到简化儒略日。
% mjd=cal2mjd(cal) 返回简化儒略日
% cal:1x6矩阵,6列分别为年月日时分秒。构造cal时可以省略末尾的0
jd = cal2jd(cal);
mjd = jd - 2400000.5;
% FUNCTION: SECONDS TO HOURS, MINUTES and SECONDS ++++++++++++++++
function timestr=secs2hms(SECS)
HOURS=SECS/60/60;
hours=floor(HOURS);
MINUTES=(HOURS-hours)*60;
minutes=floor(MINUTES);
seconds=(MINUTES-minutes)*60;
HH=num2str(hours);
MM=num2str(minutes);
SS=num2str(seconds);
if seconds<10;
SS=['0' num2str(SS)];
else
SS=num2str(SS);
end
if minutes〈10;
MM=['0’ num2str(MM)];
else
MM=num2str(MM);
end
timestr=[HH ':’ MM ':' SS];
%+++++++++++++++++++++++++++++++++++++
function timestr=cal2et(t_start,secs)
%根据开始时间,时长(秒),计算结束时间
%t_start 1×6矩阵,年,月,日,时,分,秒
%secs 时长,单位秒
t_s_days=t_start(4)/24+t_start(5)/60/24+t_start(6)/60/60/24;
t_e_days=t_s_days+secs/60/60/24;
HOURS=t_e_days*24;
hours=floor(HOURS);
展开阅读全文