gusucode.com > DSISoft是由加拿大地质调查局发布的用于垂直地震剖面(VSP)数据处理的免费软件包 > dsisoftv3/main/agcmem.m
function [dataout]=agcmem(datain,window,type) %[dataout]=agcmem(datain,window,type) % %This function will do automatic gain control with a running window equation %on the traces in datain. %The size of the sliding window is specified by the parameter 'window' in %seconds. %'type' = 1 use absolute values for normalizing %'type' = 2 use energy values (x^2) for normalizing % %Rewritten to conserve memory! %Is much slower than agc but conserves memory %Use only when memory restrictions make agc unuseable! % %Based on module by Kristen Beaty Dec. 1997 %Written by Marko Mah February 1999 %$Id: agcmem.m,v 3.0 2000/06/13 19:19:38 gilles Exp $ %$Log: agcmem.m,v $ %Revision 3.0 2000/06/13 19:19:38 gilles %Release 3 % %Revision 2.0 1999/05/21 18:45:04 mah %Release 2 % %Revision 1.1 1999/02/10 19:09:29 mah %Initial revision % % % %Copyright (C) 1998 Seismology and Electromagnetic Section/ %Continental Geosciences Division/Geological Survey of Canada % %This library is free software; you can redistribute it and/or %modify it under the terms of the GNU Library General Public %License as published by the Free Software Foundation; either %version 2 of the License, or (at your option) any later version. % %This library is distributed in the hope that it will be useful, %but WITHOUT ANY WARRANTY; without even the implied warranty of %MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %Library General Public License for more details. % %You should have received a copy of the GNU Library General Public %License along with this library; if not, write to the %Free Software Foundation, Inc., 59 Temple Place - Suite 330, %Boston, MA 02111-1307, USA. % %DSI Consortium %Continental Geosciences Division %Geological Survey of Canada %615 Booth St. %Ottawa, Ontario %K1A 0E9 % %email: dsi@cg.nrcan.gc.ca disp('[dataout]=agcmem(datain,window,type)') dataout=datain; tstart=datain.fh{9}; %start time in seconds int=datain.fh{8}; %sampling interval in seconds samples=datain.fh{7}; %number of points per trace nrec=datain.fh{12}; %number of records in datain w=round(window/int)+1; %convert 'window' from seconds to indexes pt=round(w/2); %index of point in the centre of the window % the following checks to see if the window chosen is too large if w>=samples %error check w = samples-1; pt=round(w/2); ntime = (samples - 1)*int; text = sprintf('WARNING! Window has been reset to %8.5f s',ntime); disp(text) end %if % the following checks to see if the either type 1 or 2 has been chosen if((type ~=1) & (type ~=2)) type=1; text = sprintf('WARNING! Type has been reset to 1'); disp(text) end %if % the following for loop applies the agc to each record for COUNT1=1:nrec % the following applies the method of agc specified by type switch type case 1 %absolute values % first determine how many traces there are [a,ntraces]=size(datain.dat{COUNT1}); %ntraces is the number of traces in the record trccor=zeros(samples,1); %initializes the trace correction for COUNT2=1:ntraces % first take the absolute value of the data and divide by w to make it faster temp=abs(datain.dat{COUNT1}(:,COUNT2))/w; % first determine the average in the window fact=sum(temp(1:w)); % now apply this to the first half of the window trccor(1:pt)=fact; % now apply the agc to the centre portion of the trace using a for loop for k=1:samples-w % the correction factor fact is being recalculated after each position moved fact=fact-temp(k)+temp(k+w); trccor(k+pt)=fact; end %for k % now apply this correction factor to the last half of the window i=(k+pt+1):samples; %i are the positions that still need to be corrected trccor(i)=fact; % now apply the trace correction to the data trccoreps=trccor+eps; %eps is added in to prevent divide by zero problem dataout.dat{COUNT1}(:,COUNT2)=datain.dat{COUNT1}(:,COUNT2)./trccoreps; end %for COUNT2 case 2 %squared values % first determine how many traces there are [a,ntraces]=size(datain.dat{COUNT1}); %ntraces is the number of traces in the record trccor=zeros(samples,1); %initializes the trace correction for COUNT2=1:ntraces % first square the data and divide by w squared to make it faster temp=datain.dat{COUNT1}(:,COUNT2).*datain.dat{COUNT1}(:,COUNT2)/w/w; % first determine the average in the window fact=sum(temp(1:w)); % now apply this to the first half of the window trccor(1:pt)=fact; % now apply the agc to the centre portion of the trace using a for loop for k=1:samples-w % the correction factor fact is being recalculated after each position moved fact=fact-temp(k)+temp(k+w); trccor(k+pt)=fact; end %for k % now apply this correction factor to the last half of the window i=(k+pt+1):samples; %i are the positions that still need to be corrected trccor(i)=fact; % now apply the trace correction to the data trccoreps=sqrt(trccor+eps); %eps is added in to prevent divide by zero problem dataout.dat{COUNT1}(:,COUNT2)=datain.dat{COUNT1}(:,COUNT2)./trccoreps; % the following balances the energy from trace to trace temp=dataout.dat{COUNT1}(:,COUNT2).*dataout.dat{COUNT1}(:,COUNT2); fact=sum(temp); fact=sqrt(fact); if fact==0 fact=1; %avoid divide by zero error for dead traces end %if dataout.dat{COUNT1}(:,COUNT2)=dataout.dat{COUNT1}(:,COUNT2)/fact; %applies the correction end %for COUNT2 end %type end %loop over records