conv_FD.c 2.96 KB
Newer Older
tilman.metz's avatar
tilman.metz committed
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
tilman.metz's avatar
tilman.metz committed
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 105 106 107
 *
 * 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>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Function for convolving           
 *   19/01/02, T. Bohlen
 *   last update 19.09.11 M. Schaefer => FFTW3 check http://fftw.org/
 *  ----------------------------------------------------------------------*/

#include "fd.h"
#include "segy.h"
#include <fftw3.h>

void  conv_FD(float * temp_TS, float * temp_TS1, float * temp_conv, int ns){

	/* declaration of local variables */
	int i,j, h, nfreq, npad;
	float xr, yr, x, y, dump, a;
	double npadd;
		
	/* declaration of variables for FFTW3*/
	fftw_complex *in1, *in2, *out1, *out2, *out3;
	fftw_plan p1, p2, p3;
	
    
	npad = (int)(pow(2.0, ceil(log((double)(ns))/log(2.0))+2.0) );  /* ns -> npad for usage in FFT*/
        npadd = (double)npad;
	
	in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	
	in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	
	out3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
	
	
	for (j=0;j<ns;j++){
	in1[j][0]=temp_TS[j+1];
	in1[j][1]=0.0;
		}
	for (j=ns;j<npad;j++){
	in1[j][0]=0.0;
	in1[j][1]=0.0;
		}
	
	for (j=0;j<ns;j++){
	in2[j][0]=temp_TS1[j+1];
	in2[j][1]=0.0;
		}
	for (j=ns;j<npad;j++){
	in2[j][0]=0.0;
	in2[j][1]=0.0;
		}
	
	
	p1 = fftw_plan_dft_1d(npad, in1, out1, 1, FFTW_ESTIMATE);
	fftw_execute(p1);
	
	p2 = fftw_plan_dft_1d(npad, in2, out2, 1, FFTW_ESTIMATE);
	fftw_execute(p2);
		
	/* multiplication of complex vectors */
	for(i=0;i<npad;i++){
	
	out3[i][0]    = (out1[i][0] * out2[i][0]) - (out1[i][1] * out2[i][1]);
	out3[i][1] = (out1[i][0] * out2[i][1]) + (out1[i][1] * out2[i][0]);
		   
			}
	
	p3 = fftw_plan_dft_1d(npad, out3, in1, -1, FFTW_ESTIMATE);
	fftw_execute(p3);
		
	 /* write output into data matrix */
        for(j=0;j<ns;j++){
	       temp_conv[j+1] = (1.0/npadd)*in1[j][0];
	    }
	
	
	fftw_free(in1);
	fftw_free(out1);
	
	fftw_free(in2);
	fftw_free(out2);
	
	fftw_free(out3);
	
	fftw_destroy_plan(p1);
	fftw_destroy_plan(p2);
	fftw_destroy_plan(p3);
	
}