s_spectrum.m

实现地震勘探中

switch param.plot
case 'amp'
   handle4amp_axis=gca;
   if strcmpi(param.scale,'log')
      set(handle4amp_axis,'Yscale','log')
   end
   for ii=1:nseis
      ltext{ii}=mnem2tex(datasets{ii}.name);
      plot_amplitude_spectrum_no1(s_rm_trace_nulls(datasets{ii}),param,ii);
   end
   finish_amplitude_plot_no3(param,ltext)

case 'phase'
   handle4phase_axis=gca;
   for ii=1:nseis
      plot_phase_spectrum_no2(s_rm_trace_nulls(datasets{ii}),param,ii);
      ltext{ii}=mnem2tex(datasets{ii}.name);
   end
   finish_phase_plot_no4(param,ltext)

case 'both'
   if ~strcmpi(param.figure,'new')
      error('If amplitude and phase are to be plotted the figure parameter must be set to ''new''.')
   end

   %    Create subplot axes
   if strcmpi(param.orient,'portrait')
      handle4amp_axis=subplot(2,1,1);
      handle4phase_axis=subplot(2,1,2);
   else
      handle4amp_axis=subplot(1,2,1);
      handle4phase_axis=subplot(1,2,2);
   end

   for ii=1:nseis
      ltext{ii}=mnem2tex(datasets{ii}.name);
      axes(handle4amp_axis)
      plot_amplitude_spectrum_no1(s_rm_trace_nulls(datasets{ii}),param,ii);
      axes(handle4phase_axis)
      plot_phase_spectrum_no2(s_rm_trace_nulls(datasets{ii}),param,ii);
   end
   finish_phase_plot_no4(param,ltext)
   axes(handle4amp_axis)
   finish_amplitude_plot_no3(param,ltext)

%       Link the frequency axes of the two plots
   linkaxes([handle4amp_axis,handle4phase_axis],'x')
%{
   if nargout > 0
      aux.zoom_handles=disable_zoom(figure_handle);
   else
      disable_zoom(figure_handle)
   end
%}
otherwise
   error(['Unknown plot option: "',param.plot,'".'])
end


if nargout > 0
   aux.figure_handle=figure_handle;
   try
      aux.handle4amp_axis=handle4amp_axis;
   catch  %#ok
   end
   try
      aux.handle4phase_axis=handle4phase_axis;
   catch  %#ok
   end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function aux=plot_amplitude_spectrum_no1(dataset,param,index)
% Plot amplitude spectrum

%     Apply taper window if required
switch param.window
   case {'none','rect'}
      % Do nothing
   otherwise
      dataset=s_window(dataset,param.window);
end

%       Compute FFT and amplitude spectrum
nsamp=size(dataset.traces,1);
nfft=max(nsamp,param.padding);
ftseis=fft(dataset.traces,nfft);

%       Compute frequencies values
if mod(nfft,2) == 0
   nffth=nfft/2+1;
else
   nffth=(nfft+1)/2;
end
ftseis(nffth+1:end,:)=[];
aftseis=mean(abs(ftseis),2)/nfft;

%       Normalize the spectra (if requested)
switch param.normalize
case 'peak'
   aftseis=aftseis/max(aftseis);
case 'power'
   aftseis=aftseis/norm(aftseis);
case 'amplitude'
   aftseis=aftseis/mean(aftseis);
case 'no'
   % Do nothing
otherwise
   disp(['Alert from S_SPECTRUM: Unknown normalization parameter: "',param.normalize,'". Spectrum is not normalized.'])
end


%       Set the scale to be plotted
switch lower(param.scale)
case {'power','log'}
   aftseis=aftseis.^2;
otherwise
   % do nothing
end


%       Set the scale to be plotted
%aftseis(1)=0;
switch lower(param.scale)
case 'log'
   aftseis=max(aftseis,eps*max(aftseis));
case 'db'
   aftseis=max(aftseis,eps*max(aftseis));
   aftseis=20*log10(aftseis);
   aftseis=aftseis-max(aftseis);
otherwise
   % do nothing
end

%       Determine abscissa values
nyquist=500/dataset.step;
freq=(0:2:nfft)*nyquist/nfft;
bool=freq >= param.frequencies{1} &  freq <= param.frequencies{2};
freq=freq(bool);
aftseis=aftseis(bool);

%       Plot amplitude spectrum
idx=mod(index-1,length(param.linestyles))+1;
line(freq,aftseis,'LineWidth',param.linewidth,'Color',param.colors{idx}, ...
   'LineStyle',param.linestyles{idx});

if nargout > 0
   aux.ftseis=ftseis;
   aux.freq=freq;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function aux=plot_phase_spectrum_no2(dataset,param,index)
% Plot phase spectrum

[nsamp,ntr]=size(dataset.traces);
nfft=max(param.padding,nsamp);

switch param.timezero
case 'best'
   timezero=[];
case 'actual'
   timezero=dataset.first;
otherwise
   error(['Unknown value of parameter "timezero": ',param.timezero])
end

if ntr > 1
  for ii=1:ntr
      [phase,aux]=signal_phase_spectrum(dataset.traces(:,ii),nfft, ...
                  timezero,param.tiepoint);
      if ii==1;
         avphase=phase;
      else
         avphase=avphase+phase;
      end
   end
   phase=avphase/ntr;
else
   [phase,aux]=signal_phase_spectrum(dataset.traces,nfft, ...
        timezero,param.tiepoint);
end

%       Determine abscissa values
nyquist=500/dataset.step;
freq=aux.freq*nyquist;
%freq=(0:2:nfft)*nyquist/nfft;
bool=freq >= param.frequencies{1} &  freq <= param.frequencies{2};
freq=freq(bool);
phase=phase(bool);
ik=true;
while all(phase < 0) && any(phase < 180)
   phase=phase+360;
   ik=false;
end
if ik
   while all(phase > 0) && any(phase > 180)
      phase=phase-360;
   end
end

%       Plot phase spectrum
idx=mod(index-1,length(param.linestyles))+1;
line(freq,phase,'LineWidth',param.linewidth,'Color',param.colors{idx}, ...
   'LineStyle',param.linestyles{idx});

if nargout > 0
   aux.phase=phase;
   aux.freq=freq;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function finish_amplitude_plot_no3(param,ltext)

%       Create axis labels
xlabel(param.label4x)

switch lower(param.scale)
case 'linear'
   switch lower(param.normalize)
   case {'peak','power'}
      label4y='Normalized spectral amplitude';
   otherwise
      label4y='Spectral amplitude';
   end

case 'power'
   switch lower(param.normalize)
   case {'peak','power'}
      label4y='Normalized spectral power';
   otherwise
      label4y='Spectral power';
   end

case 'log'
   switch lower(param.normalize)
   case {'peak','power'}
      label4y='Normalized spectral power';
   otherwise
      label4y='Spectral power';
   end

case 'db'
   switch lower(param.normalize)
   case {'peak','power'}
      label4y='Normalized spectral power (dB)';
   otherwise
      label4y='Spectral power (dB)';
   end

otherwise
   label4y='';
end

%{
switch param.normalize
case 'peak'
   label4y=[label4y,' (normalized peak amplitude)'];
case 'power'
   label4y=[label4y,' (normalized energy)'];
otherwise
   % do nothing
end
%}

ylabel(label4y)

grid on
box on
bgGray
axis tight

legend(ltext,'Location',param.lloc)

initiate_2d_tracking({'Freq.',param.xunits,'Frequency'},{'Amp.','n/a','Amplitude'})


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function finish_phase_plot_no4(param,ltext)

%       Create axis labels
xlabel(param.label4x)

if isempty(param.timezero)
   ylabel('Phase in degrees (zero-time at peak of Hilbert transform)')
else
   ylabel('Phase in degrees')
end

grid on
box on
bgGray
axis tight

legend(ltext,'Location',param.lloc)

initiate_2d_tracking({'Freq.',param.xunits,'Frequency'},{'Phase.','degrees','Phase'})