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
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