Classical Gaussian filter
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=R8), | intent(in), | dimension(1:long, 1:larg) | :: | tab |
2D array in |
|
| integer(kind=I4), | intent(in) | :: | long |
2D array width |
||
| integer(kind=I4), | intent(in) | :: | larg |
2D array height |
||
| real(kind=R8), | intent(in) | :: | cutoff |
cut-off wavelength |
||
| real(kind=R8), | intent(out), | dimension(1:long, 1:larg) | :: | bf_tab |
2D array out |
|
| logical(kind=I4), | intent(in) | :: | multi_fft |
multiple fft at once ? |
||
| real(kind=R8), | intent(in), | optional | :: | pad |
fft padding |
|
| character(len=*), | intent(in), | optional | :: | ext |
extension |
|
| character(len=*), | intent(in), | optional | :: | type_apo |
apodization type |
subroutine fft_filter(tab, long, larg, cutoff, bf_tab, multi_fft, pad, ext, type_apo) !================================================================================================ !! Classical Gaussian filter !------------------------------------------------------------------------------------------------ use surfile, only : init_scal, write_surf, SCALE_SURF implicit none integer(kind=I4), intent(in ) :: long !! *2D array width* integer(kind=I4), intent(in ) :: larg !! *2D array height* real (kind=R8), intent(in ) :: cutoff !! *cut-off wavelength* logical(kind=I4), intent(in ) :: multi_fft !! *multiple fft at once ?* real (kind=R8), intent(in ), optional :: pad !! *fft padding* character(len=*), intent(in ), optional :: ext !! *extension* character(len=*), intent(in ), optional :: type_apo !! *apodization type* real (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab !! *2D array in* real (kind=R8), intent(out), dimension(1:long, 1:larg) :: bf_tab !! *2D array out* integer(kind=I4) :: nx2, ny2, iex, iey, ibx, iby real (kind=R8) :: o_pad logical(kind=I4) :: with_pad character(len=:), allocatable :: o_ext character(len=:), allocatable :: o_type_apo complex(kind=R8), dimension(:,:), allocatable :: cmple real (kind=R8), dimension(:,:), allocatable :: tab_ext, gauss_tab, tab_tmp with_pad = .true. if ( .not.present( ext ) ) then o_ext = 'symmetry' else o_ext = ext endif if ( .not.present( type_apo ) ) then o_type_apo = 'no_apo' else o_type_apo = type_apo endif if ( .not.present( pad ) ) then o_pad = PAD_FFT_FILTER else if ( pad < 0. ) then o_pad = 1 nx2 = long ny2 = larg with_pad = .false. else o_pad = pad endif endif if ( with_pad ) then nx2 = 2 * ( nint(o_pad * long)/2 ) ny2 = 2 * ( nint(o_pad * larg)/2 ) endif allocate( tab_ext(1:nx2, 1:ny2) ) ! allocate( tab_tmp(1:nx2, 1:ny2) ) ! allocate( cmple(1:nx2, 1:ny2) ) ! if ( nx2 > long ) then ibx = ceiling( (nx2 - long)/2. ) ; iex = ibx + long - 1 iby = ceiling( (ny2 - larg)/2. ) ; iey = iby + larg - 1 call extend( tab_in = tab(1:long, 1:larg), & ! tab_out = tab_ext(1:nx2, 1:ny2), & ! nx = long, & ! ny = larg, & ! nx2 = nx2, & ! ny2 = ny2, & ! ext = o_ext, & ! type_apo = o_type_apo ) ! else tab_ext(1:nx2, 1:ny2) = tab(1:long, 1:larg) endif if (multi_fft) then call tab_calc_fftw3_real_fwd( tab_in = tab_ext(1:nx2, 1:ny2), & ! tab_ou = cmple(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! else call calc_fftw3_real_fwd( tab_in = tab_ext(1:nx2, 1:ny2), & ! tab_ou = cmple(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! endif allocate( gauss_tab(1:nx2, 1:ny2) ) call gaussian_filter( long = nx2, & ! larg = ny2, & ! xc = cutoff, & ! gauss_filt = gauss_tab(1:nx2, 1:ny2) ) ! cmple(1:nx2, 1:ny2) = cmple(1:nx2, 1:ny2) * gauss_tab(1:nx2, 1:ny2) deallocate( gauss_tab ) if (multi_fft) then call tab_calc_fftw3_real_bwd( tab_in = cmple(1:nx2, 1:ny2), & ! tab_ou = tab_ext(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! else call calc_fftw3_real_bwd( tab_in = cmple(1:nx2, 1:ny2), & ! tab_ou = tab_ext(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! endif if ( nx2 > long ) then bf_tab(1:long, 1:larg) = tab_ext(ibx:iex, iby:iey) else bf_tab(1:long, 1:larg) = tab_ext(1:nx2, 1:ny2) endif deallocate(cmple) deallocate(tab_ext, tab_tmp) return endsubroutine fft_filter