s_spectrum.m

基于Matlab的地震数据处理显示和测井数据显示于处理的小程序

         atext='Power';
      end
      amin=0;

    else
      atext='Amplitude (dB)';
      amp=amp/max(amp);
%      idx=find(amp < 1.0e-5);     
      amp(amp < 1.0e-5)=1.0e-5;   % Constrain the possible values of the amplitude spectrum
      amp=20*log10(amp/max(amp));
      amin=-inf;
    end

  end

  if phaseplot
             % Account for start time of signal 
     phshift=exp(0.002*pi*i*f*stemp.first)';
     if ntr > 1
       temp=mean(ft(1:nffth,:).*phshift(:,ones(ntr,1)),2);
     else
       temp=ft(1:nffth,:).*phshift;
     end  
 
     phase=atan2(imag(temp),real(temp));
     if ntr > 0 
       if strcmpi(param.option,'average')  
          phase=mean(phase,2);
       elseif strcmpi(param.option,'envelope')
          alert(' There is no meaningful interpretation of option "envelope" for the phase')
          alert(' average used instead')
          phase=mean(phase,2);
       else
         disp([' Unknown option ',param.option])
         disp(' Passible values are: "envelope" and "average"')
         error(' Abnormal termination')
       end
     end
     if unwrap_phase
        phase=unwrap(phase);         % Unwrap phase (MATLAB Signal Processing Toolbox)
     end
     phase=phase*(180/pi);
  end

  if param.frequencies{1} > 0 | param.frequencies{2} < endfreq;
%    idx=find(f > param.frequencies{1} & f < param.frequencies{2});
    ff=[param.frequencies{1},f(f > param.frequencies{1} & ...
                               f < param.frequencies{2}),param.frequencies{2}];
    if ampplot
       amp=interp1(f,amp,ff,'*linear');
       if strcmpi(param.scale,'linear') & strcmpi(param.normalize,'yes')
          amp=amp/max(amp);
       end
    end
    if phaseplot
       phase=interp1(f,phase,ff,'*linear');
    end
    f=ff;
  end

%               Create legend
  if ntr == 1
    attrib='';
  else
    attrib=attribute;
  end
  if ~isempty(dsetname)
    ltext(ii)={[strrep(dsetname,'_','\_'),attrib]};
  else
    ltext(ii)={['Input data ',num2str(ii),attrib]};
  end

  if ampplot		% Scale amplitude spectrum if desired
    if (strcmpi(param.scale,'linear') | strcmpi(param.scale,'power')) & strcmpi(param.normalize,'yes')
       amp=amp/max(amp);
    end
    if strcmpi(param.scale,'power')
       amp=amp.^2;
    end
  end

  if ampplot & phaseplot
     hh1=subplot(2,1,1); 
     bgGray
     amp_handles(ii)=plot(f,amp,get_color(ii,param.colors),'LineWidth',param.linewidth);
     if ii == 1
        axis([param.frequencies{1},param.frequencies{2},amin,inf])
     end
     hold on
     subplot(2,1,2); 
     phase_handles(ii)=plot(f,phase,get_color(ii,param.colors),'LineWidth',param.linewidth);
     bgGray
     hold on
     if unwrap_phase
        if nseis == 1
           ll=floor(min(phase/100))*100;
           uu=ceil(max(phase/100))*100;
           axis([param.frequencies{1},param.frequencies{2},ll,uu]);
        end

     elseif ii == 1
        axis([param.frequencies{1},param.frequencies{2},-180,180])
     end

  elseif ampplot
     amp_handles(ii)=plot(f,amp,get_color(ii,param.colors),'LineWidth',param.linewidth);
     bgGray
     hold on

  elseif phaseplot
     phase_handles(ii)=plot(f,phase,get_color(ii,param.colors),'LineWidth',param.linewidth);
     bgGray
     hold on
  else
    error(' Neither amplitude spectrum nor phase spectrum selected')
  end   
     
end

if strcmp(param.figure,'new')
   timeStamp;
end

units=getfield(varargin{1},'units');

if ampplot & phaseplot

  if isempty(param.lloc)     % Set location of legend
     loc=5;
  else
     loc=param.lloc;
  end
  if isempty(param.legend)
     legend(char(ltext),loc);
  else
     legend(param.legend,loc)
  end
  title('Phase Spectrum')
  if strcmpi(units,'ms')
     xlabel('Hz')
  elseif strcmpi(units,'s')
     xlabel('mHz')
  elseif strcmpi(units,'m')
     xlabel('Wavelengths per 1000 m')
  elseif strcmpi(units,'ft')
     xlabel('Wavelengths per 1000 ft')
  elseif strcmpi(units,'samples')
     xlabel('Wavelengths per 1000 samples')
  else
    % Unknown units: do nothing
  end
  ylabel('Phase angle in degree')
  grid on, zoom on
 
  axes(hh1)
  if isempty(param.legend)
     legend(char(ltext),loc);
  else
     legend(param.legend,loc)
  end
  
  title('Amplitude Spectrum')
  if strcmpi(units,'ms')
     xlabel('Hz')
  elseif strcmpi(units,'s')
     xlabel('mHz')
  elseif strcmpi(units,'m')
     xlabel('Wavelengths per 1000 m')
  elseif strcmpi(units,'ft')
     xlabel('Wavelengths per 1000 ft')
  elseif strcmpi(units,'samples')
     xlabel('Wavelengths per 1000 samples')
  else
    % Unknown units: do nothing
  end
  ylabel(atext)
  grid on, zoom on

elseif ampplot
  if isempty(param.lloc)     % Set location of legend
     loc=1;
  else
     loc=param.lloc;
  end
  if isempty(param.legend)
     legend(char(ltext),loc);
  else
     legend(param.legend,loc)
  end
  title('Amplitude Spectrum')

  if strcmpi(units,'ms')
     xlabel('Hz')
  elseif strcmpi(units,'s')
     xlabel('mHz')
  elseif strcmpi(units,'m')
     xlabel('Wavelengths per 1000 m')
  elseif strcmpi(units,'ft')
     xlabel('Wavelengths per 1000 ft')
  elseif strcmpi(units,'samples')
     xlabel('Wavelengths per 1000 samples')
  else
    % Unknown units: do nothing
  end
  ylabel(atext)
  grid on, zoom on

else
 if isempty(param.lloc)     % Set location of legend
     loc=5;
  else
     loc=param.lloc;
  end
  if isempty(param.legend)
     legend(char(ltext),loc);
  else
     legend(param.legend,loc)
  end
  title('Phase spectrum')

  if strcmpi(units,'ms')
     xlabel('Hz')
  elseif strcmpi(units,'s')
     xlabel('mHz')
  elseif strcmpi(units,'m')
     xlabel('Wavelengths per 1000 m')
  elseif strcmpi(units,'ft')
     xlabel('Wavelengths per 1000 ft')
  elseif strcmpi(units,'samples')
     xlabel('Wavelengths per 1000 samples')
  else
    % Unknown units: do nothing
  end
  
  ylabel('Phase angle in degree')
  grid on, zoom on

end
hold off

if ampplot
   aux.amp_handles=amp_handles;
end
if phaseplot
   aux.phase_handles=amp_handles;
end

if ampplot & phaseplot
   linkedzoom('onx2');
   if nargout > 0
      aux.zoom_handles=disable_zoom(aux.figure_handle);
   else
      disable_zoom(gcf)
   end
end


if nargout == 0
   clear aux
else
   if ampplot
      aux.amp_handles=amp_handles;
   end
   if phaseplot
      aux.phase_handles=amp_handles;
   end
end