timedomain_filt.c 3.6 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
Florian Wittkamp's avatar
Florian Wittkamp committed
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
Tilman Steinweg's avatar
Tilman Steinweg committed
5
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
Tilman Steinweg's avatar
Tilman Steinweg committed
7
8
9
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
Tilman Steinweg's avatar
Tilman Steinweg committed
11
12
13
14
15
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
Florian Wittkamp's avatar
Florian Wittkamp committed
16
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
Tilman Steinweg's avatar
Tilman Steinweg committed
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Filter seismograms in time domain with a Butterworth filter
 *   Lowpass or highpass filtering can be applied                                  
 *  ----------------------------------------------------------------------*/
#include "fd.h"
#include "segy.h"
#include "cseife.h"

void  timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int method){

	/* 
	data	: 	2-dimensional array containing seismograms (
	fc	:	corner frequency in Hz
	order	:	order of filter
	ntr	:	number of traces
	ns	:	number of samples
	method	:	definition of filter
			1: lowpass filter
			2: highpass filter
	*/

	/* declaration of extern variables */
	extern float DT, F_HP;
	extern int ZERO_PHASE, NT,MYID;
	
	/* declaration of local variables */
	int itr, j, ns_reverse;
	double *seismogram, *seismogram_reverse, T0;
47
	double *seismogram_hp, *seismogram_reverse_hp, T0_hp;
Tilman Steinweg's avatar
Tilman Steinweg committed
48
49
50
51
	
	seismogram = dvector(1,ns);
	if (ZERO_PHASE==1) seismogram_reverse = dvector(1,ns);
	
52
53
54
	seismogram_hp = dvector(1,ns);
	if (ZERO_PHASE==1) seismogram_reverse_hp = dvector(1,ns);
	
Tilman Steinweg's avatar
Tilman Steinweg committed
55
	T0=1.0/(double)fc;
56
57
58
59
60
61
62
63
64
65
66
	if(F_HP)
		T0_hp=1.0/(double)F_HP;
	if(method==2)
		T0_hp=1.0/(double)fc;
		
	if (method==1){    /* lowpass filter */
		for (itr=1;itr<=ntr;itr++){
			for (j=1;j<=ns;j++){
				seismogram[j]=(double)data[itr][j];
			}
			
67
			seife_lpb(seismogram,ns+1,DT,T0,order); /* ns+1 because vector[0] is also allocated and otherwise seife_lpb do not filter the last sample */
68
			
Tilman Steinweg's avatar
Tilman Steinweg committed
69
70
71
72
73
74
75
76
77
78
79
80
			if (ZERO_PHASE==1){
			ns_reverse=ns;
				for (j=1;j<=ns;j++) {
					seismogram_reverse[ns_reverse]=seismogram[j];
					ns_reverse--;}
			seife_lpb(seismogram_reverse,ns+1,DT,T0,order);
			ns_reverse=ns; 
				for (j=1;j<=ns;j++) {
					seismogram[ns_reverse]=seismogram_reverse[j];
					ns_reverse--;}
			}

81
82
83
84
85
			for (j=1;j<=ns;j++){
				data[itr][j]=(float)seismogram[j];
			}
		}
	} /* end of itr<=ntr loop */
Tilman Steinweg's avatar
Tilman Steinweg committed
86

87
88
89
90
91
	if ((method==2)||(F_HP)){   /*highpass filter*/
		for (itr=1;itr<=ntr;itr++){
			for (j=1;j<=ns;j++){
				seismogram_hp[j]=(double)data[itr][j];
			}
Tilman Steinweg's avatar
Tilman Steinweg committed
92
			
93
			seife_hpb(seismogram_hp,ns+1,DT,T0_hp,order);
Tilman Steinweg's avatar
Tilman Steinweg committed
94
95
96
97
			
			if (ZERO_PHASE==1){
			ns_reverse=ns;
				for (j=1;j<=ns;j++) {
98
					seismogram_reverse_hp[ns_reverse]=seismogram_hp[j];
Tilman Steinweg's avatar
Tilman Steinweg committed
99
					ns_reverse--;}
100
			seife_hpb(seismogram_reverse_hp,ns+1,DT,T0_hp,order);
Tilman Steinweg's avatar
Tilman Steinweg committed
101
102
			ns_reverse=ns;
				for (j=1;j<=ns;j++) {
103
					seismogram_hp[ns_reverse]=seismogram_reverse_hp[j];
Tilman Steinweg's avatar
Tilman Steinweg committed
104
105
					ns_reverse--;}
			}
106
107
108
			for (j=1;j<=ns;j++){
				data[itr][j]=(float)seismogram_hp[j];
			}
Tilman Steinweg's avatar
Tilman Steinweg committed
109
110
111
112
113
114
		}
	} /* end of itr<=ntr loop */
	
	free_dvector(seismogram,1,ns);
	if (ZERO_PHASE==1) free_dvector(seismogram_reverse,1,ns);

115
116
117
	free_dvector(seismogram_hp,1,ns);
	if (ZERO_PHASE==1) free_dvector(seismogram_reverse_hp,1,ns);

Tilman Steinweg's avatar
Tilman Steinweg committed
118
} /* end of function */