收藏 分销(赏)

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

上传人:二*** 文档编号:4536390 上传时间:2024-09-27 格式:DOC 页数:67 大小:1.30MB 下载积分:5 金币
下载 相关 举报
基于Matlab的TEQC绘图程序代码doc.doc_第1页
第1页 / 共67页
本文档共67页,全文阅读请下载到手机保存,查看更方便
资源描述
基于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);
展开阅读全文

开通  VIP会员、SVIP会员  优惠大
下载10份以上建议开通VIP会员
下载20份以上建议开通SVIP会员


开通VIP      成为共赢上传

当前位置:首页 > 学术论文 > 其他

移动网页_全站_页脚广告1

关于我们      便捷服务       自信AI       AI导航        抽奖活动

©2010-2025 宁波自信网络信息技术有限公司  版权所有

客服电话:0574-28810668  投诉电话:18658249818

gongan.png浙公网安备33021202000488号   

icp.png浙ICP备2021020529号-1  |  浙B2-20240490  

关注我们 :微信公众号    抖音    微博    LOFTER 

客服