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: Handel 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.compiled % 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=logical(1);
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*zeros(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
