s_frequency_attributes.m

实现地震勘探中

function attributes=s_frequency_attributes(seismic,varargin)
% Compute median frequency, bandwidth, and power spectrum of seismic input data
% Written by: E. R.: February 22, 2005
% Last updated: March 7, 2005: Handle also matrices
%
%              attributes=s_frequency_attributes(seismic,varargin)
% INPUT
% seismic      seismic data or matrix
% varargin     one or more cell arrays; the first element of each cell array is a keyword,
%              the other elements are parameters. Presently, keywords are:
%     'window'   window to use for spectrum computation. Possible window names are:
%	      'Hamming', 'Hanning', 'Nuttall',  'Papoulis', 'Harris',
%	      'Rect',    'Triang',  'Bartlett', 'BartHann', 'Blackman'
%	      'Gauss',   'Parzen',  'Kaiser',   'Dolph',    'Hanna'.
%	      'Nutbess', 'spline',  'none'
%             'Rect' and 'none' have the same effect as has {'window',[]};
%             Case is irrelevant (e.g. 'Hanning' and 'hanning' are the same)
%             If the number of traces goes into the thousands then windowing 
%             becomes time consuming and may not be necessary.
%             Default: {'window','none'}
% OUTPUT
% attributes   structure with fields
%              'bandwidth'         average bandwidth of the seismic data set
%              'bandwidth_octaves' bandwidth in octaves
%              'median_frequency'  median frequency
%              'power_spectrum'    power spectrum
%              'frequencies'       frequency associated with each sample 
%                                  of "power_spectrum"
%             If the input data set is a matrix, frequencies and bandwidth 
%             are in fractions of the Nyquist frequency
%             If no output argument is supplied the function prints out 
%             bandwidth and median frequency.

global ABORTED S4M

param.window=[];
param.average='yes';

param=assign_input(param,varargin);

if isnumeric(seismic)  &&  size(seismic,1) > 1 % "seismic" is matrix}
   seismic=s_convert(seismic,0,500);
end
    
[nsamp,ntr]=size(seismic.traces);
nfft=max(2^nextpow2(nsamp),256);

%     Compute bandwidth, median frequency, and integrated spectrum
if isempty(param.window)  ||  strcmpi(param.window,'none')  ||  strcmpi(param.window,'rect')
   ftraces=fft(seismic.traces,nfft);
else
%       Set up display of progress bar
   gui_active(1);       % Add an abort button
   if ~S4M.deployed     % Reset persistent variable in "progressbar"
      clear progressbar % Commented out to allow compilation
   end

   h=progressbar([],0,'Computing windowed spectrum ...');
   wind=mywindow(nsamp,param.window);
   ftraces=zeros(nfft,ntr);
   for ii=1:ntr
      if ishandle(h)
         h=progressbar(h,1/ntr);
      end
      if ~gui_active
         ABORTED=true;
         progressbar(h,-1);
         return
      end
      drawnow

      ftraces(:,ii)=fft(seismic.traces(:,ii).*wind,nfft);
   end
   progressbar(h,-1);
end

nfreq=nfft/2+1; 	                     % Length(ftraces) is even by design
power_spectrum=abs(ftraces(1:nfreq,:)).^2;   % Keep only positive frequencies 

if strcmpi(param.average,'yes')
   power_spectrum=sum(power_spectrum,2);
   ntr=1;
end

cs=cumsum(power_spectrum);

fnyquist=500/seismic.step;
df=2*fnyquist/nfft;
freq=(0:1:nfreq-1)*df;


med_freq=NaN(1,ntr);
for ii=1:ntr
   if cs(end,ii) > 0
      index=find(diff(cs(:,ii)) > 0);
      med_freq(ii)=interp1(cs(index,ii),freq(index),cs(end,ii)*0.5);
   end
end

%	Compute bandwidth
bandwidth=NaN*zeros(1,ntr);
temp=max(power_spectrum);
idx=find(temp > 0);
bandwidth(idx)=df*cs(end,idx)./temp(idx);  % Mean spectrum divided by spectral maximum
octaves=log2((med_freq+0.5*bandwidth)/(max(med_freq-0.5*bandwidth,eps)));

if nargout > 0
   attributes.bandwidth=bandwidth;
   attributes.bandwidth_octaves=octaves;
   attributes.median_frequency=med_freq;
   attributes.power_spectrum=power_spectrum;
   attributes.frequencies=freq;
else
   disp(['Bandwidth = ',num2str(round(bandwidth)),' Hz'])
   disp(['Bandwidth in octaves = ',num2str(round(100*octaves)/100)])
   disp(['Median frequency = ',num2str(round(med_freq)),' Hz'])
end

%figure
%plot(freq,power_spectrum)
%keyboard