Hi,

I am trying to use rlft3 to filter a 2D image. I just want to apply a Gaussian filter after applying an FFT and then inv FFT to recover the filtered image. I have been struggling for the implementation, although I expected it be quite straightforward. I would be really really grateful if someone could have a quick look at the short fortran code below and shed me some light. Note I fft the gaussian filter, although it is not necessary. I was just trying to understand the structure a bit better.

************************
program filter

implicit none

!uses rlft3
integer n1,n2,n3,i,j
parameter (n1=32,n2=32,n3=1)
real data(n1,n2),dx,dy,pi
real fil(n1,n2),xx,yy,sigma,lambda,xl,yl,pixl,piyl,fi,f j
complex spec(n1/2,n2),speq(n2),spec2(n1/2,n2),speq2(n2)
equivalence (data,spec)
equivalence (fil,spec2)


speq=0.
pi=atan(1.)*4.
sigma=1.
dx=.25
dy=.25
xl=(n1-1)*dx
yl=(n2-1)*dy
lambda=xl/2.

open(7,file='data1.txt',status='unknown')

! a topography should loaded
do j=1,n2
do i=1,n1
call random_number(xx)
xx=0.
data(i,j)=(xx*2.-1)+sin(2.*pi*float(i-1)*dx/lambda) !+sin(2.*pi*float(j-1)*dy/lambda)
fil(i,j)=1./sqrt(2.*pi)/sigma*exp(-((float(i-1)*dx)**2+(float(i-1)*dy)**2)/2./sigma**2)
write(7,*) float(i-1)*dx,float(j-1)*dy,data(i,j)
enddo
enddo
fil=fil/sum(fil)
close(7)

! computes the real FFT in 2D
print*,'computing the fft'
call rlft3(data,speq,n1,n2,n3,1)
call rlft3(fil,speq2,n1,n2,n3,1)

!applies filter
print*,'applying filter'
do j=1,n2
do i=1,n1
!real part
data(i,j)=data(i,j)*fil(i,j)-data(i+1,j)*fil(i+1,j)
!imaginary part
! data(i+1,j)=data(i+1,j)*fil(i,j)+fil(i+1,j)*data(i ,j)
enddo
enddo
!spec=spec*spec2

! computes the inv FFT to load the filtered topography
print*,'computing the ifft'
call rlft3(data,speq,n1,n2,n3,-1)
data=data/float(n1*n2*n3)*2.

open(8,file='data2.txt',status='unknown')
do j=1,n2
do i=1,n1
write(8,*) float(i-1)*dx,float(j-1)*dy,data(i,j)
enddo
enddo
close(8)

Thanks,
fred