convolution using fft without equivalence


mohan
10-01-2008, 10:52 AM
Hello,

I modified the segment of example code for convolving two functions using FFT, since I have problems with equivalence statement for various reasons. However, the result doesnt make sense. I am doing something silly, I am sure, probably in unwrapping the result of rlft3 ... any help would be appreciated !

Thanks
Mohan

image1 and image2 are nXm arrays.

call rlft3(image1,speq1,n,m,1,-1)
call rlft3(image2,speq2,n,m,1,-1)

image1=image1*2/(n*m)
image2=image2*2/(n*m)

do j=1,m
speq(j)=2.d0/(n*m)*speq1(j)*speq2(j)
end do


do i=1,n,2
do j=1,m
indy=j
indx=(i-1)/2+1
realim1(indx,indy)=image1(i,j)
imagim1(indx,indy)=image1(i+1,j)
realim2(indx,indy)=image2(i,j)
realc(indx,indy)=realim1(indx,indy)*realim2(indx,i ndy)-
imagim1(indx,indy)*imagim2(indx,indy)
imagc(indx,indy)=realim1(indx,indy)*imagim2(indx,i ndy)+
realim2(indx,indy)*imagim1(indx,indy)

conv(i,j)=realc(indx,indy)
conv(i+1,j)=imagc(indx,indy)
end do
end do

call rlft3(conv,speq,n,m,1,1)

do i=1,n,2
do j=1,m
indy=j
indx=(i-1)/2+1
realc(indx,indy)=conv(i,j)
imagc(indx,indy)=conv(i+1,j)
end do
end do

mohan
10-05-2008, 09:16 AM
Hello,

Before any of you reply, I did manage to fix the problem and I get the correct answer now.

However, the accuracy of rlft3 seems to be about 1 in 10^8 whereas fftw has an accuracy of 1 in 10^16 or so. Is that normal ?

Mohan