harmon.m,v

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

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


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

2.0
date	99.05.21.18.45.45;	author mah;	state Exp;
branches;
next	1.2;

1.2
date	99.01.29.16.29.52;	author adam;	state Exp;
branches;
next	1.1;

1.1
date	99.01.06.19.09.03;	author kay;	state Exp;
branches;
next	;


desc
@@


3.0
log
@Release 3
@
text
@function [dataout] = harmon(datain,freq1,freq2,inter)
%  function [dataout] = harmon(datain,freq1,freq2,inter)
%
%  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.
%
%  Usage:
%
%  	Input:
%		datain  : Input data
%		freq1   : Start scanning at this frequency (Hz)
%		freq2   : Upper frequency limit (Hz)
%		inter   : Frequency step (Hz)
%
%	Output:
%		dataout : Output data
%
%	Example :
%		This is what I use to remove 60 Hz noise
%		filtered = harmon(data,59.88,60.12,0.01);
%
%	Reference:
%	         http://www.cg.NRCan.gc.ca/staff/adam/software/monofreq.html
%
%written by:  Erick ADAM
%  Dec. 1995

%$Id: harmon.m,v 2.0 1999/05/21 18:45:45 mah Exp gilles $
%$Log: harmon.m,v $
%Revision 2.0  1999/05/21 18:45:45  mah
%Release 2
%
%Revision 1.2  1999/01/29 16:29:52  adam
%Updated to work with Dsisoft variables and add comments
%
%Revision 1.1  1999/01/06 19:09:03  kay
%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(datain,freq1,freq2,inter)')
deltat=datain.fh{8};      %sampling interval in seconds
freq = freq1:inter:freq2; %Vector of frequencies for DFT
freqs = length(freq);     %Number of frequencies for DFT
samples = datain.fh{7};   %number of points per trace;

%Initialize local variables for speedup
cos_wt = ones(samples,freqs);   % Cosine lookup table
sin_wt = ones(samples,freqs);   % 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 : freqs,
	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

% 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
 		large= 0.0; % Maximum amplitude (reset to 0 before each scan)
		for j = 1: freqs,  % Loop over frequencies
			trace = datain.dat{COUNT}(:,tr);
			suma = sum(trace .* cos_wt(:,j))* 2.0 /samples;
			sumb = sum(trace .* sin_wt(:,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
                % 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
@d29 1
a29 1
%$Id: harmon.m,v 1.2 1999/01/29 16:29:52 adam Exp mah $
d31 3
@


1.2
log
@Updated to work with Dsisoft variables and add comments
@
text
@d29 1
a29 1
%$Id: harmon.m,v 1.1 1999/01/06 19:09:03 kay Exp $
d31 3
@


1.1
log
@Initial revision
@
text
@d1 2
a2 2
function [filtered] = harmon(data,freq1,freq2,inter,deltat)
%  function [filtered] = harmon(data,freq1,freq2,inter,deltat)
d11 4
a14 5
%		data :  Matrix of data to filter
%		freq1 : Start scanning at this frequency
%		freq2 : Upper frequency limit
%		inter : Frequency step (Hz)
%		deltat : Sample rate
d17 1
a17 1
%		Filtered : Matrix of filtered data
d21 1
a21 1
%		filtered = harmon(data,59.88,60.12,0.01,0.0001);
d23 3
d29 4
a32 2
%$Id:$
%$Log:$
d34 1
d62 5
a66 1
disp('[filtered] = harmon(data,freq1,freq2,inter,deltat)')
d68 6
a73 8
freq = freq1:inter:freq2;
freqs = length(freq);
samples = size(data,1);
traces = size(data,2);
cos_wt = ones(samples,freqs);
sin_wt = ones(samples,freqs);
filtered = ones(samples,traces);
samp = (0.0 : deltat : (samples-1)*deltat)';
d76 3
a78 3
	w = 6.283185307 * freq(i);
	cos_wt(:,i)=cos(samp*w);
	sin_wt(:,i)=sin(samp*w);
a79 18
% Noise estimation
for i = 1 :traces,
	large= 0.0;
	for j = 1: freqs,
		suma = sum(data(:,i) .* cos_wt(:,j))* 2.0 /samples;
		sumb = sum(data(:,i) .* sin_wt(:,j))* 2.0 /samples;
		amp_max = (suma * suma) + (sumb*sumb);
		if(amp_max > large)
			freqn=freq(j);
			freqi=j;
			large = amp_max;
			a_max=suma;
			b_max=sumb;
		end
	end
% The filter is applied here
	filtered(:,i) = data(:,i)-((a_max.*cos_wt(:,freqi))+(b_max .* sin_wt(:,freqi)));
end
d81 20
d102 1
@