filespcgrm.m

2001年美国南中国海试验数据文件.专门用来读数据的.

% ---------------------------------------------------------------------
%		FILESPCGRM
%  Script for plotting a spectrogram from a KVDH datafile.
%  Can pick an individual channel, or choose to average over several
%  channels.  Currently set to sample once every datarecord, which for
%  a typical file gives 300 samples over a 2hour period, about once every
%  30 seconds.
%
% -----------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% routine: filespcgrm.m
% author: B. Sperry
% date: 6/10/98
% history:
%
%		hacked on hard by Newhall 
%
% notes:
%   o Files needed:  getAllDRH.m, getDRH.m, getSeqs.m, wfall.m, db20.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Needed in case files are shortened
clear T data F P fdata


% cd /home/tomo/arthur/Primer3/shark/noise
% cd  /home/tomo/arthur/AsiaEx01/prelim


% --- Parameters (to set as you wish...) ---
%
 
% AsiaEx files
filename = './05150249.CSD';
filename = './05150109.CSD';

  chns     = 3	;		% channel numbers (will average if more than 1)
  slen     = 1024;	% this is also the FFT-size
  nseqs    = 5;		% No. of FFT sequences to average
  Tstart   = 0.0;	% time to start sampling datefile [minutes]
  Tend     = 10;	% time to halt sampling (-1 for end of file) [minutes]
  Tsamp    = -1;	% interval at which to sample datafile [minutes]
                        % Tsamp = -1 to sample once every data record
  % Plot parameters
  plot_type = 0;        % 0 - color spec'gram, 1 - waterfall plot
  dbmin     = -50;	% dB levels for plotting
  dbmax     = 50;
  scaling   = '';       % select 'auto' if want automatic dB scaling

  % these are optional, they will be set automatically. Use
  % only when you really want the plots to stay a certain way!!                      
  YMAX = 800;           % Max Frequency used in plotting
  YMIN = 0;             % min freq.
  
%-----------------------------------------------
% Shouldn't need to change anything below here
% ----------------------------------------------


% --- Open file and get DRH information. ---

  fid   = fopen(filename,'rb','ieee-le');
  if(fid<0), error(['Error opening datafile : ' filename]); end;
  drhs  = getAllDRH(fid);
  nrecs = length(drhs);
  fsamp = drhs(1).rhfs;

% --- Set up time vector for sampling datafile ---

  t_begin = drhs(1).time(1)   + drhs(1).time(2)/60000;
  t_end   = drhs(end).time(1) + drhs(end).time(2)/60000 + ...
            (drhs(end).date(2)-drhs(1).date(2))*24*60 + ...
	    drhs(end).npts/drhs(end).rhfs/60; 
  tlen_datafile = t_end-t_begin;
  if(Tend<=0), Tend=tlen; end;
  if(Tsamp<=0),
    sampvec = [0:-1:-(nrecs-1)];
  else,
    sampvec = [Tstart:Tsamp:Tend]*60.0;
  end;
  nsamps = length(sampvec);

% --- Read seqs in, FFT'ing as we go... ---

  fprintf('Reading from datafile: %s\n',filename);
  %% fprintf('Internal datafile name: %s\n',drhs(1).fname);
  fdata = zeros(slen/2+1,nsamps);
  win   = repmat(hanning(slen),[1,nseqs]);
  F     = [0:slen/2]/slen*fsamp;
  
  save_data = [];
  for k=1:nsamps,
    [data,t,gt]=getSeqs(fid, sampvec(k),  nseqs, slen);
    P = fft(squeeze(mean(data(chns,:,:),1)).*win,[],1);
    fdata(:,k) = mean(abs(P(1:slen/2+1,:)),2);
    T(k)=gt(1).min+gt(1).msec/60000;
    fprintf('%g\n',T(k)/60);
  end;
  fclose(fid);
  clear P data win 

% --- Plotting ---

  if(strcmp(scaling,'auto')),
    dbmax = round(max(db20(fdata(:))));
    dbmin = dbmax - 50;
  end;
  
  
  if(plot_type==0),
    fig=figure('Name','SPECGRM');
    clf;
    ax(1)=gca;
    imagesc(T/60, F, db20(fdata), [dbmin dbmax] );
    grid on;
    set(gca,'YDir','normal'); 
    colormap(jet);
    
    % reference lines???
    %% h=hline([136 145 256 265 336 345 448 457 556 565], ...
    %% 	'linestyle',':','color','k','LineWidth',1.0);

    xlabel('Time (hrs)','FontWeight','Bold'); 
    ylabel('Frequency (Hz)','FontWeight','Bold'); 
    if (exist('YMAX')),
	set(gca,'Ylim',[YMIN,YMAX]);
    end
    
    tstr=sprintf('%s %s',drhs(1).aexp,drhs(1).proj);
    title(tstr,'FontWeight','Bold');
    hc=colorbar; 
    set(get(hc,'Title'),'string','dB');
    set(hc,'Position',[0.83 0.11 0.029 0.4075]);
    set(hc,'Ylim',[dbmin dbmax]);

    
    % --- Text Axes ---
    
    ax(1)=gca;
    ax(2)=axes('Position',[0.80 0.60 0.03 0.4],'Visible','off');
    ht(1)=text(0.1,.75,drhs(1).adate);
    ht(2)=text(0.1,.65,drhs(1).atime);
    ht(3)=text(0.1,.55,drhs(1).fname);
    ht(4)=text(0.1,.4,sprintf('Fs: %g Hz',drhs(1).rhfs));
    ht(5)=text(0.1,.3,sprintf('FFT size: %d',slen));
    ht(6)=text(0.1,.2,'Window: hanning');
    ht(7)=text(0.1,.1,sprintf('FFT avgs: %d',nseqs));
    s=sprintf('%1d ',chns-1);
    ht(8)=text(0.1,0,[ 'Chns: ' s ]);
    
    set(ht(1:3),'FontName','Times');
    set(ht(4:8),'FontName','Times','FontSize',9);
    
    axes(ax(1));
    orient('portrait');

	% k = find(filename == '.');
	% filename2 = filename(1:k-1);
	% S = sprintf('print -depsc -epsi %s.eps',filename2)
	% eval(S);

	% disp('   Paused....  hit <return>');
	% pause
	% disp('            OK... continuing ');

	%% axis([ 8 8.5 100 500]);
	%% S = sprintf('print -depsc -epsi %s_in.eps',filename2)
	%% eval(S);



  elseif(plot_type==1),    % WATERFALL plotting

    % some optional indexing for plotting a subsection
    %indf = find(200<=F & F<=500);
    %indt = find(0<=T & T<=18.5*60);
    indt = [1:length(T)];
    indf = [1:length(F)];
    
    wfall(F(indf),db20(fdata(indf,indt)).',T(indt)/60,[dbmin dbmax], ...
	'fill','Color','b','gain',8,'ydir','normal');

    orient('tall'); 
    xlabel('Freq (Hz)'); ylabel('Time (hrs)');
    tstr=sprintf('%s %s: %s, %s',drhs(1).aexp,drhs(1).proj, ...
	drhs(1).adate,drhs(1).atime);
    title(tstr,'FontWeight','Bold');
    orient('tall');
    
  else,
  
    fprintf('Sorry, we don''t do those kinds of plots around here.');
    
  end;

% ---------------------------------------------------------------------