| real*8 trig(20),x(0:9)
real*8 work(10)
n=10
m=1
do i=1,n
x(i-1)=2.d0*i
end do
ifail=0
do j=1,n
write(*,*)'input x,index ', x(j-1),j
end do
c
CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
do j=1,n
write(*,*)'output x,index ', x(j-1),j
end do
end
OUTPUT
---------
rawhide.rto.dec.com> a.out
input x,index 2.00000000000000 1
input x,index 4.00000000000000 2
input x,index 6.00000000000000 3
input x,index 8.00000000000000 4
input x,index 10.0000000000000 5
input x,index 12.0000000000000 6
input x,index 14.0000000000000 7
input x,index 16.0000000000000 8
input x,index 18.0000000000000 9
input x,index 20.0000000000000 10
output x,index 2.00000000000000 1
output x,index 20.0000000000000 2
output x,index 18.0000000000000 3
output x,index 16.0000000000000 4
output x,index 14.0000000000000 5
output x,index 12.0000000000000 6
output x,index 10.0000000000000 7
output x,index 8.00000000000000 8
output x,index 6.00000000000000 9
output x,index 4.00000000000000 10
|
|
Well, at least something is correct. Data is scaled by one of the two
routines.
The output from C06FQF has not been properly bit reversed. Please see whether
the pair C06EAF and C06EBF correctly recover the original sequence. If they
do file a QAR against C06FQF. If they don't, it's likely that these routine
cannot do forward and then reverse in the sense we know. I suspect repeating
the call pair an extra time will produce the original data:
eg. do
CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
instead of
CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
If the little gimmick above works, it means the trig tables for C06FPF and
C06FQF are the same, when one should be the conjugate of the other.
Meanwhile, what does the program do with the data after the FFT?
|
|
>The output from C06FQF has not been properly bit reversed. Please see whether
>the pair C06EAF and C06EBF correctly recover the original sequence.
Unfortunately I haven't got them
>I suspect repeating
>the call pair an extra time will produce the original data:
>eg. do
> CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
> call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
> CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
> call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
>
>instead of
>
> CALL C06FPF(m,n,x,'I',TRIG,WORK,IFAIL)
> call c06fqf(m,n,x,'I',TRIG,WORK,IFAIL)
Absolutely right
>Meanwhile, what does the program do with the data after the FFT?
SUBROUTINE FOURIER(F,FT,NN1,NN2,TRIG1,TRIG2,WORK,INIT)
C calculates the 2-d Fourier transform of
C array F(B)=Psi(A+B)*Psi(A-B) using their Symmetry
C
C Link with NAGDP/LIB.
IMPLICIT NONE
CHARACTER*1 INIT
INTEGER N1,N2,N(2),NH(2),NN1,NN2,I,K,IHILF,KHILF,IFAIL
REAL*8 F(0:NN1/2*NN2-1),FT(0:NN2*NN1-1),
& TRIG1(2*NN1),TRIG2(2*NN2),WORK(NN1*NN2)
EXTERNAL C06FPF,C06FQF
N(1)=NN1
N(2)=NN2
NH(1)=N(1)/2
NH(2)=N(2)/2
C initialize to nil:
DO 1 I=0,N(2)-1
DO 2 K=0,N(1)-1
C FT(I,K)=0.D0
FT(K*N(2)+I)=0.D0
2 CONTINUE
1 CONTINUE
C first FFT (rows ). F is real therefore these are Hermit transforms
C rows
IFAIL=0
CALL C06FPF(NH(1),N(2),F(0),INIT,TRIG2,WORK,IFAIL)
C change order of the array elements for the second FFT (columns )
C (siehe Text S. 16):
I=0
KHILF=0
DO 3 K=0,NH(1)-1
IHILF=K
C FT(I,K)=F(IHILF,KHILF)
FT(K*N(2)+I)=F(KHILF*NH(1)+IHILF)
3 CONTINUE
DO 4 I=1,NH(2)-1
KHILF=I
DO 5 K=0,NH(1)-1
IHILF=K
C FT(I,K)=F(IHILF,KHILF)
FT(K*N(2)+I)=F(KHILF*NH(1)+IHILF)
C if(abs(F(KHILF*NH(1)+IHILF)) .gt. 1.0D-5) then
C write(6,*)F(KHILF*NH(1)+IHILF)
C endif
5 CONTINUE
KHILF=N(2)-I
DO 6 K=NH(1)+1,N(1)-1
IHILF=N(1)-K
C FT(I,K)=F(IHILF,KHILF)
FT(K*N(2)+I)=F(KHILF*NH(1)+IHILF)
6 CONTINUE
4 CONTINUE
I=NH(2)
KHILF=I
DO 7 K=0,NH(1)-1
IHILF=K
C FT(I,K)=F(IHILF,KHILF)
FT(K*N(2)+I)=F(KHILF*NH(1)+IHILF)
7 CONTINUE
DO 8 I=NH(2)+1,N(2)-1
KHILF=N(2)-I
DO 9 K=0,NH(1)-1
IHILF=K
C FT(I,K)=F(IHILF,KHILF)
FT(K*N(2)+I)=F(KHILF*NH(1)+IHILF)
9 CONTINUE
KHILF=I
DO 10 K=NH(1)+1,N(1)-1
IHILF=N(1)-K
C FT(I,K)=-F(IHILF,KHILF)
FT(K*N(2)+I)=-F(KHILF*NH(1)+IHILF)
C if(abs(F(KHILF*NH(1)+IHILF)) .gt. 1.0D-13) then
C write(6,*)abs(F(KHILF*NH(1)+IHILF))
C endif
10 CONTINUE
8 CONTINUE
C second FFT (columns=rows of FT). FT is hermit, therefore
C the transformed rows are real:
CALL C06FQF(N(2),N(1),FT(0),INIT,TRIG1,WORK,IFAIL)
END
|