calc_envelope.c 3.22 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
/*-----------------------------------------------------------------------------------------
 * Copyright (C) 2013  For the list of authors, see file AUTHORS.
 *
 * This file is part of DENISE.
 * 
 * DENISE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
 * DENISE 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 General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Calculate the envelope of a seismogram via Hilbert transform                                  
 *   last update 27/1/10, L. Rehor
 *   update to FFTW3 01/07/2013 Martin Schaefer
 *  ----------------------------------------------------------------------*/
#include "fd.h"
#include <fftw3.h>

void calc_envelope(float ** datatrace, float ** envelope, int ns, int ntr){

	/* declaration of variables */
	int i,j, nfreq, npad, k;
	float xr, yr, x, y, dump, a, *h;
	double npadd;
		
	/* declaration of variables for FFTW3*/
	fftw_complex *in, *out;
	fftw_plan p1, p2;

	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
	/* calculation of the Hilbert transform */
	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */


	/* calculation of the FFT */
	npad = (int)(pow(2.0, ceil(log((double)(ns))/log(2.0))+2.0) );  /* ns -> npad for usage in FFT*/
        npadd = (double)npad;
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	
	/* calculation of the vector h */
	     h = vector(1,npad);
	     if ((2*(npad/2))==npad) {
		/* npad is even */
		h[1]=1.0;
		h[npad/2+1]=1.0;
		for (i=2;i<=(npad/2);i++) {
	   	     h[i] = 2.0; }
	     } else {
		/* npad is odd */
		h[1]=1.0;
		for (i=2;i<=((npad+1)/2);i++) {
	   	     h[i]=2.0;}
	     }
	     
	for(k=1;k<=ntr;k++){
						
	     for (j=0;j<ns;j++){
	     in[j][0]=datatrace[k][j+1];
	     in[j][1]=0.0;
		     }
	     for (j=ns;j<npad;j++){
	     in[j][0]=0.0;
	     in[j][1]=0.0;
		     }
	     /* FFT */
	     p1 = fftw_plan_dft_1d(npad, in, out, 1, FFTW_ESTIMATE);
	     fftw_execute(p1);
	     fftw_destroy_plan(p1);		    

	     /* elementwise multiplication of FFT and h */
	     for (i=0;i<npad;i++){
		out[i][0] *= h[i+1];
		out[i][1] *= h[i+1];
	     }


	     /* inverse FFT */
	     p2 = fftw_plan_dft_1d(npad, out, in, -1, FFTW_ESTIMATE);
	     fftw_execute(p2);
             fftw_destroy_plan(p2);

	     /* %%%%%%%%%%%%%%%%%%%%%%% */
	     /* calculation of envelope */
	     /* %%%%%%%%%%%%%%%%%%%%%%% */
	     for (j=0;j<ns;j++){
	     	envelope[k][j+1]=(1.0/npadd)*sqrt(in[j][0]*in[j][0]+in[j][1]*in[j][1]);
	     }
	     
	}

	fftw_free(in);
	fftw_free(out);
	free_vector(h,1,npad);
	
}