基于Matlab的TEQC绘图程序代码

1、%在 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 电离层延迟观测值变化率(米/秒)

2、 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

3、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);%按行读取文件至数组AA=importdata(path file,&

4、#39;t');%定义SAT，存放卫星数据%GPS有32颗卫星，存放序号1-32，%GLONASS有32颗卫星，存放序号33-64，%BEIDOU有35颗卫星，存放序号65-99SAT(1:length(A),1:99)=NaN;%sats(1:length(A),1:99)=NaN;%存放采样时间，单位秒tsec(1:length(A)=NaN;%读取文件的第一行filelx=A1;%判断是哪种格式switch filelx case 'COMPACT' %读取数据采样间隔 t_samp=char(A(3); %读取开始时间 mjl=char(A(4); %读取数据采

5、样间隔 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=

6、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_

7、start = deblank(t_start); s = splitstr(t_start,'*',6); %t_start_time=char(s2) '年' char(s3) '月' char(s4) '日' char(s5) '时' char(s6) '分' num2str(str2num(char(s7),'%02d') '秒' ; %获取采样的开始时间，2013,12，7,03,05,55 t_s_time=str2num(char(s2), str2nu

8、m(char(s3),str2num(char(s4),str2num(char(s5),str2num(char(s6),str2num(char(s7); n=1;i=3; otherwise disp('数据格式存在问题'); returnend%n=1;i=3;%sats=str2num(Ai);snyggfilen=strrep(filen,'_','-');%生成进度条h=waitbar(0,'正在读取数据，请稍等 ' char(snyggfilen);while i<length(A); waitbar(i/l

9、ength(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(s1); if s2='-

10、1' %如果为-1，使用上次的卫星编号行 satbhstr=oldsatbh; s = splitstr(satbhstr,'*',32); else oldsatbh=satbhstr; end %获取当前卫星数量 satcount=str2num(s2); case 'COMPACT','COMPACT2' %记录采样时间，单位秒 tsec(n)=(n-1)*T_SAMP; if s1='-1' %如果为-1，使用上次的卫星编号行 satbhstr=oldsatbh; s = splitstr(satbhstr,

11、9;*',32); else oldsatbh=satbhstr; end %获取当前卫星数量 satcount=str2num(s1); 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(sk-1); case 'COMPACT

12、2' %如果是COMPACT2格式，卫星编号行第2个数据开始是卫星编号 satbh=char(sk-1); case 'COMPACT3' %如果是COMPACT3格式，卫星编号行第3个数据开始是卫星编号 satbh=char(sk); end switch satbh(1); case 'G' %如果是GPS卫星，获得与编号对应的存储序号，范围1-32 index=str2num(satbh(2:3); case 'R' %如果是GLONASS卫星，获得与编号对应的存储序号，范围33-64 index=32+str2num(satbh(

13、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(sdatak-2); %如果卫星数据不是数值型，用0替代 if isempty(str2num(sdatastr)=1 sdatastr='0.000' end %将卫星数据存入SAT数组的对应位置 SAT(n,i

14、ndex)=str2num(sdatastr); end %for n=n+1; %根据数据文件结构，一行是卫星编号，下一行就是对应的卫星数据 %程序一次读取2行数据 i=i+2;end %while%数据读取完毕，关闭进度条close(h);%获取已使用的卫星编号，返回给sat_bhsat_bh=getsatbh(SAT);%增加最后一个结束字符sat_bh_temp=sat_bh' 'End'sat_bh_end=sat_bh_temp'%删去所有都为NaN值的列，返回给sat_datasat_data=delnancol(SAT);%删去所有NaN值的采样

15、时间列t_seconds=delnancol(tsec);%删去所有都为NaN值的行，返回sat_datassat_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 numberswitch filelx case 'COMPACT','COMPACT2' t_s_jd

16、_time=MJD_START; case 'COMPACT3' t_s_jd_time=datenum(t_s_time);endfor 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 onrow,col=size(sat_datas);%绘制渐变彩色图pcolor(t_jd_time

17、9;,1:col,sat_datas');set(gcf,'renderer','zbuffer');shading flatset(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);s

18、et(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(

19、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);%+%以

20、下是本函数中所使用的一些相关函数%+%+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 stringindex=temp; ind

21、ex=index+1; temp=; end else temp=temp,str(i); end endelse %常规分隔符 for i=1:length(str) if str(i)=split if isempty(temp)=false stringindex=temp; index=index+1; temp=; end else temp=temp,str(i); end endendif isempty(temp)=false stringindex=temp;ends=string;%+function B=delnancol(A);%删除矩阵中所有值都是nan的列x y=s

22、ize(A);for i=y:-1:1 if all(isnan(A(:,i)=1 A(:,i)=; endendB=A;%+function B=delnanrow(A);%删除矩阵中所有值都是nan的行x y=size(A);for i=x:-1:1 if all(isnan(A(i,:)=1 A(i,:)=; endendB=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)=

23、0; end endendB=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',nu

24、m2str(i-32,'%02d'); end end sindex=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

25、 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 '

26、;.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&#

27、39; 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='卫星方位角(

28、6;)' 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

29、 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('vert

30、ical'); 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('horizonta

31、l'); 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：1

32、x6矩阵，6列分别为年月日时分秒。构造cal时可以省略末尾的0% 公元1582年10月4日24:00点之前使用儒略历，公元1582年10月15日00:00点之后使用公历if length(cal) < 6cal(6)=0;endyear=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 tha

33、t there is no year 0 A.D.endm=month;if( m <= 2 )% January and February come after December.m = m+12; y = y - 1;ende=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=15

34、82& month=10 &day<15)b = -38;elseb = floor(a/4) - a); % number of century years that are not leap yearsendc=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

35、日24:00点之前使用儒略历，公元1582年10月15日00:00点之后使用公历if (J < 1721423.5)% 公元1月1日0时BC = 1;elseBC = 0;end%start from Julian March 1, 4801 B.C.if( J < 2299160.5 )% before 1582.10.4. 24:00 is Julian calenderj0=floor(J+0.5);dd=J+0.5-j0;else% after 1582.10.15. 00:00 is Gregorian calender%number of certury years t

36、hat are not leap yearj0=n1+n2+n3+J+10;dd=j0+0.5-floor(j0+0.5);j0=floor(j0+0.5);endj0=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>12month=month-12;year=year+1;endyear=year-BC;sec=round(dd*86400);hour

37、=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时可以省略末尾的0jd = 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-minute