KJ135
/
ISCE_INSAR
mirror of http://172.16.0.12:5000/WBJY/ISCE_INSAR.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

use FFTW to speed up filtering

LT1AB
CunrenLiang 2020-04-24 21:42:04 -07:00 committed by piyushrpt
parent e96a110221
commit 9794a75f55
1 changed files with 36 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
contrib/alos2filter/src/psfilt1.c
View File

@ -2,6 +2,7 @@
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include <fftw3.h>
#define PLUS 1 #define PLUS 1
#define MINU 2 #define MINU 2
@ -75,6 +76,11 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw
FILE *int_file, *sm_file; FILE *int_file, *sm_file;
double psd,psd_sc; /* power spectrum, scale factor */
int ii, jj;
fftwf_plan p_forward;
fftwf_plan p_backward;
int k; int k;
float sf; // scale factor for FFT, otherwise FFT will magnify the data by FFT length float sf; // scale factor for FFT, otherwise FFT will magnify the data by FFT length
// usually the magnitude of the interferogram is very large, so the data are // usually the magnitude of the interferogram is very large, so the data are
@ -212,6 +218,9 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw
} }
lc=0; lc=0;
p_forward = fftwf_plan_dft_2d(fftw, fftw, (fftw_complex *)seg_fft[0], (fftw_complex *)seg_fft[0], FFTW_FORWARD, FFTW_MEASURE);
p_backward = fftwf_plan_dft_2d(fftw, fftw, (fftw_complex *)seg_fft[0], (fftw_complex *)seg_fft[0], FFTW_BACKWARD, FFTW_MEASURE);
for (i=0; i < yh; i += step){ for (i=0; i < yh; i += step){
for(i1=fftw - step; i1 < fftw; i1++){ for(i1=fftw - step; i1 < fftw; i1++){
fread((char *)cmp[i1], sizeof(fcomplex), width, int_file); fread((char *)cmp[i1], sizeof(fcomplex), width, int_file);
@ -229,8 +238,30 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw
if(i%(2*step) == 0)fprintf(stderr,"\rprogress: %3d%%", (int)(i*100/yh + 0.5)); if(i%(2*step) == 0)fprintf(stderr,"\rprogress: %3d%%", (int)(i*100/yh + 0.5));
for (j=0; j < width; j += step){ for (j=0; j < width; j += step){
psd_wgt(cmp, seg_fft, alpha, j, i, fftw); //psd_wgt(cmp, seg_fft, alpha, j, i, fftw);
fourn((float *)seg_fft[0]-1,nfft,ndim,isign); /* 2D inverse FFT of region, get back filtered fringes */ ////////////////////////////////////////////////////////////////////////////////////////
//replace function psd_wgt with the following to call FFTW, crl, 23-APR-2020
for (ii=0; ii < fftw; ii++){ /* load up data array */
for (jj=j; jj < j+fftw; jj++){
seg_fft[ii][jj-j].re = cmp[ii][jj].re;
seg_fft[ii][jj-j].im = cmp[ii][jj].im;
}
}
//fourn((float *)seg_fft[0]-1, nfft, ndim, -1); /* 2D forward FFT of region */
fftwf_execute(p_forward);
for (ii=0; ii < fftw; ii++){
for (jj=0; jj < fftw; jj++){
psd = seg_fft[ii][jj].re * seg_fft[ii][jj].re + seg_fft[ii][jj].im * seg_fft[ii][jj].im;
psd_sc = pow(psd,alpha/2.);
seg_fft[ii][jj].re *= psd_sc;
seg_fft[ii][jj].im *= psd_sc;
}
}
/////////////////////////////////////////////////////////////////////////////////////////
//fourn((float *)seg_fft[0]-1,nfft,ndim,isign); /* 2D inverse FFT of region, get back filtered fringes */
fftwf_execute(p_backward);
for (i1=0; i1 < fftw; i1++){ /* save filtered output values */ for (i1=0; i1 < fftw; i1++){ /* save filtered output values */
for (j1=0; j1 < fftw; j1++){ for (j1=0; j1 < fftw; j1++){
@ -285,6 +316,9 @@ int psfilt1(char *inputfile, char *outputfile, int width, double alpha, int fftw
fprintf(stdout,"\nnumber of lines written to file: %d\n",lc); fprintf(stdout,"\nnumber of lines written to file: %d\n",lc);
stop_timing(); stop_timing();
fftwf_destroy_plan(p_forward);
fftwf_destroy_plan(p_backward);
free(bufcz); free(bufcz);
//free_cmatrix(cmp, 0, fftw-1, -fftw,width+fftw); //free_cmatrix(cmp, 0, fftw-1, -fftw,width+fftw);
//free_cmatrix(sm, 0,fftw-1,-fftw,width+fftw); //free_cmatrix(sm, 0,fftw-1,-fftw,width+fftw);