harmon_new.m,v

具有特色的地震数据处理源码

head	3.0;
access;
symbols;
locks; strict;
comment	@// @;


3.0
date	2000.06.13.19.20.30;	author gilles;	state Exp;
branches;
next	2.0;

2.0
date	99.05.21.18.45.46;	author mah;	state Exp;
branches;
next	1.10;

1.10
date	99.02.10.20.47.00;	author mah;	state Exp;
branches;
next	1.9;

1.9
date	99.02.10.20.11.36;	author perron;	state Exp;
branches;
next	1.8;

1.8
date	99.02.10.20.07.19;	author perron;	state Exp;
branches;
next	1.7;

1.7
date	99.02.10.20.04.11;	author perron;	state Exp;
branches;
next	1.6;

1.6
date	99.02.10.19.52.34;	author perron;	state Exp;
branches;
next	1.5;

1.5
date	99.02.10.16.15.19;	author eadam;	state Exp;
branches;
next	1.4;

1.4
date	99.02.10.13.50.49;	author mah;	state Exp;
branches;
next	1.3;

1.3
date	99.02.02.15.55.12;	author mah;	state Exp;
branches;
next	1.2;

1.2
date	99.02.02.15.07.07;	author mah;	state Exp;
branches;
next	1.1;

1.1
date	99.02.02.14.20.47;	author mah;	state Exp;
branches;
next	;


desc
@modified harmon to use a specified time window to calculate the corrections
@


3.0
log
@Release 3
@
text
@function [dataout] = harmon_new(datain,freq1,freq2,inter,t1,t2)
%  function [dataout] = harmon_new(datain,freq1,freq2,inter,t1,t2)
%
%  Harmon is a time domain notch filter.  Discrete Fourier
%  Transforms are calculated to estimate noise phase and amplitude.
%  The estimated noise is subtracted from the trace in the time domain.
%  It has been rewritten to analyze a window of data from t1 to t2
%  and then to apply it to the whole trace
%
%  Usage:
%
%  	Input:
%		datain  : Input data
%		freq1   : Start scanning at this frequency (Hz)
%		freq2   : Upper frequency limit (Hz)
%		inter   : Frequency step (Hz) such that ((freq2-freq1)/inter) is
%                          a multiple of 10
%               t1      : Start time of window (seconds)
%               t2      : End time of window (seconds)
%
%	Output:
%		dataout : Output data
%
%	Example :
%		This is what I use to remove 60 Hz noise
%		filtered = harmon_new(data,59.88,60.12,0.01,0,0.5);
%
%	Reference:
%	         http://www.cg.NRCan.gc.ca/staff/adam/software/monofreq.html
%
% Warning: Program does not check for reasonable parameters
%
%written by:  Erick ADAM
%Dec. 1995
%Rewritten by Marko Mah February 1999


%$Id: harmon_new.m,v 2.0 1999/05/21 18:45:46 mah Exp gilles $
%$Log: harmon_new.m,v $
%Revision 2.0  1999/05/21 18:45:46  mah
%Release 2
%
%Revision 1.10  1999/02/10 20:47:00  mah
%fixing up bugs previous programmer(s) put in
%
%Revision 1.9  1999/02/10 20:11:36  perron
%*** empty log message ***
%
%Revision 1.8  1999/02/10 20:07:19  perron
%modifying if statement at line 119
%
%Revision 1.7  1999/02/10 20:04:11  perron
%Modifying if statement at line 119
%
%Revision 1.6  1999/02/10 19:52:34  perron
%Modifying help text
%
%Revision 1.5  1999/02/10 16:15:19  eadam
%Optimisation by iterative windowing fro finding the maximum
%
%Revision 1.4  1999/02/10 13:50:49  mah
%fixed bug
%
%Revision 1.3  1999/02/02 15:55:12  mah
%corrected amplitude of correction
%
%Revision 1.2  1999/02/02 15:07:07  mah
%fixed a bug
%
%Revision 1.1  1999/02/02 14:20:47  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] = harmon_new(datain,freq1,freq2,inter,t1,t2)')
deltat=datain.fh{8};      %sampling interval in seconds
freq = freq1:inter:freq2; %Vector of frequencies for DFT
freqlen = length(freq);     %Number of frequencies for DFT 
samples = datain.fh{7};   %number of points per trace;
% the following converts t1 and t2 to points
t1pt=round(t1/deltat)+1;
t2pt=round(t2/deltat)+1;

%Initialize local variables for speedup
cos_wt = ones(samples,freqlen);   % Cosine lookup table
sin_wt = ones(samples,freqlen);   % Sine lookup table
dataout = datain;               % Output variable
samp = (0.0 : deltat : (samples-1)*deltat)'; % Time vector

% Optimization : I create a lookup table for Sine and cosine.
for i = 1 : freqlen,
	w = 6.283185307 * freq(i); % Convert frequency to radian (w = 2*PI*f)
	cos_wt(:,i)=cos(samp*w);   % Calculate cos(w*t) and store in lookup table
	sin_wt(:,i)=sin(samp*w);   % Calculate sin(w*t) and store in lookup table
end

%step = 10; 

% Calculate amplitudes from DFT at each frequencies of vector freq
for COUNT=1:datain.fh{12} % Loop over over number of records
 for tr=1:datain.th{COUNT}(12,1)  % loop over traces
  step = floor(freqlen/10);   % Speedup to calculate DFT for every 10th frequencies 
  large= 0.0; % Maximum amplitude (reset to 0 before each scan)
  trace = datain.dat{COUNT}(t1pt:t2pt,tr);
  start_freq = 1;         % start at the first frequency
  end_freq = freqlen;     % End at the last frequency 

  if step >= 1 	  % loop until every frequency is read

   for j = start_freq:step:end_freq,  % Loop over frequency range 
    suma = sum(trace.*cos_wt(t1pt:t2pt,j))* 2.0 /samples;
    sumb = sum(trace.*sin_wt(t1pt:t2pt,j))* 2.0 /samples;
    amp_max = (suma*suma)+(sumb*sumb); %Amplitude for freq(j)

    if amp_max > large %if Amplitude freq(j) is larger than previous ones
     freqi=j; %record the index of the maximum amplitude  
     large = amp_max; % reset large
     a_max=suma;  %record factor A of the DFT
     b_max=sumb;  %record factor B of the DFT
    end % End if statement

   end % loop over frequency

   start_freq = (freqi-2)*step; % New frequency range as a function of the maximum amplitude
   end_freq = (freqi-2)*step;
   step = floor(step/10);  % reduce the step by a factor of 10
  end % End if statement
 % the following corrects a_max and b_max because a smaller window was used to calculate them
 a_max=a_max*samples/(t2pt-t1pt+1);
 b_max=b_max*samples/(t2pt-t1pt+1);
 % The filter is applied here
 dataout.dat{COUNT}(:,tr) = datain.dat{COUNT}(:,tr)-((a_max*cos_wt(:,freqi))+(b_max * sin_wt(:,freqi)));
 end % loop over traces
end % loop over records


@


2.0
log
@Release 2
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.10 1999/02/10 20:47:00 mah Exp mah $
d40 3
@


1.10
log
@fixing up bugs previous programmer(s) put in
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.9 1999/02/10 20:11:36 perron Exp $
d40 3
@


1.9
log
@*** empty log message ***
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.8 1999/02/10 20:07:19 perron Exp perron $
d40 3
d118 1
a118 1
step = 10; 
d123 1
a123 1
  step = floor(freqlen/10);   % Speedup : calculate DFT for every 10th frequencies 
d128 4
a131 2
  if (step => 1) 	  % loop until every frequency is read
   for j = start_freq: step:end_freq,  % Loop over frequency range 
d135 2
a136 1
    if(amp_max > large) %if Amplitude freq(j) is larger than previous ones
d142 1
d144 1
@


1.8
log
@modifying if statement at line 119
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.7 1999/02/10 20:04:11 perron Exp perron $
d40 3
d125 1
a125 1
  if (step => 1) % loop until every frequency is read
@


1.7
log
@Modifying if statement at line 119
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.6 1999/02/10 19:52:34 perron Exp perron $
d40 3
d122 1
a122 1
  if step >= 1 % loop until every frequency is read
@


1.6
log
@Modifying help text
@
text
@d38 1
a38 1
%$Id: harmon_new.m,v 1.5 1999/02/10 16:15:19 eadam Exp perron $
d40 3
d119 1
a119 1
  if step => 1 % loop until every frequency is read
@


1.5
log
@Optimisation by iterative windowing fro finding the maximum
@
text
@d17 1
a17 1
                          a multiple of 10
d38 1
a38 1
%$Id: harmon_new.m,v 1.4 1999/02/10 13:50:49 mah Exp eadam $
d40 3
@


1.4
log
@fixed bug
@
text
@d16 2
a17 1
%		inter   : Frequency step (Hz)
d38 1
a38 1
%$Id: harmon_new.m,v 1.3 1999/02/02 15:55:12 mah Exp mah $
d40 3
d84 1
a84 1
freqlen = length(freq);     %Number of frequencies for DFT
d103 2
d108 1
d111 18
a128 11
  for j = 1: freqlen,  % Loop over frequencies
   suma = sum(trace.*cos_wt(t1pt:t2pt,j))* 2.0 /samples;
   sumb = sum(trace.*sin_wt(t1pt:t2pt,j))* 2.0 /samples;
   amp_max = (suma*suma)+(sumb*sumb); %Amplitude for freq(j)
   if(amp_max > large) %if Amplitude freq(j) is larger than previous ones
    freqi=j; %record the index of the maximum amplitude  
    large = amp_max; % reset large
    a_max=suma;  %record factor A of the DFT
    b_max=sumb;  %record factor B of the DFT
   end % End if statement
  end % loop over frequency
@


1.3
log
@corrected amplitude of correction
@
text
@d37 1
a37 1
%$Id: harmon_new.m,v 1.2 1999/02/02 15:07:07 mah Exp mah $
d39 3
d80 1
a80 1
freqs = length(freq);     %Number of frequencies for DFT
d87 2
a88 2
cos_wt = ones(samples,freqs);   % Cosine lookup table
sin_wt = ones(samples,freqs);   % Sine lookup table
d93 1
a93 1
for i = 1 : freqs,
d104 6
a109 6
  for j = 1: freqs,  % Loop over frequencies
   suma = sum(trace.* cos_wt(t1pt:t2pt,j))* 2.0 /samples;
   sumb = sum(trace.* sin_wt(t1pt:t2pt,j))* 2.0 /samples;
   amp_max = (suma * suma) + (sumb * sumb); % Amplitude for freq(j)
   if(amp_max > large) % if Amplitude freq(j) is larger than previous ones
    freqi=j;   % record the index of the maximum amplitude  
d111 2
a112 2
    a_max=suma;  % record factor A of the DFT
    b_max=sumb;  % record factor B of the DFT
@


1.2
log
@fixed a bug
@
text
@d37 1
a37 1
%$Id: harmon_new.m,v 1.1 1999/02/02 14:20:47 mah Exp mah $
d39 3
d112 3
@


1.1
log
@Initial revision
@
text
@d30 2
d37 4
a40 2
%$Id: Exp $
%$Log:  $
d43 1
d78 1
a78 1
t2pt=rount(t2/deltat)+1;
@