This is the code for a FFT in Fortran 77. As you can see I used the subroutines provided by the chapter 12 of the book. I tried to make the fourier trasform of the function exp(-t^2) for t =(-1,1), trasform which is in terms of frecquency is √(pi)*exp{-(1/4)f²},≥0. Below the code I put the output of the code.

It something wrong with my program , the example? Do I use corect the subroutines?

Code:

program tfourier
parameter (n=16,isign=1)
real data(n),k,l
integer i,n
do i=1,n
k=i
l=(-1+2*(k-1)/(n-1))*(-1+2*(k-1)/(n-1))
data(i)=exp(-l)
enddo
call realft(data,n,isign)
do i=1,n
write(*,*) data(i)
enddo
end program

subroutine realft(data,n,isign)
integer isign,n
real data(n)
integer i,i1,i2,i3,i4,n2p3
real c1,c2,h1i,h1r,h2i,h2r,wis,wrs
double precision theta,wi,wpi,wpr,
* wr,wtemp
theta=3.141592653589793d0/(dble(n/2))
c1=0.5
if (isign.eq.1) then
c2=-0.5
call four1(data,n/2,1)
else
c2=0.5
theta=-theta
endif
wpr=-2.0d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=n+3
do i=2,n/4
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=sngl(wr)
wis=sngl(wi)
h1r=c1*(data(i1)+data(i3))
h1i=c1*(data(i2)-data(i4))
h2r=-c2*(data(i2)+data(i4))
h2i=c2*(data(i1)-data(i3))
data(i1)=h1r+wrs*h2r-wis*h2i
data(i2)=h1i+wrs*h2i+wis*h2r
data(i3)=h1r-wrs*h2r+wis*h2i
data(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+ wtemp*wpi+wi
enddo
if (isign.eq.1) then
h1r=data(1)
data(1)=h1r+data(2)
data(2)=h1r-data(2)
else
h1r=data(1)
data(1)=c1*(h1r+data(2))
data(2)=c1*(h1r-data(2))
call four1(data,n/2,-1)
endif
return
end

SUBROUTINE four1(data,nn,isign)
INTEGER isign,nn
REAL data(2*nn)
INTEGER i,istep,j,m,mmax,n
REAL tempi,tempr
DOUBLE PRECISION wr,theta,wi,wpi,wpr,wtemp
n=2*nn
j=1
do i=1,n,2
if (j.gt.i) then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=nn
1 if((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
goto 1
endif
j=j+m
enddo
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta=6.28318530717959d0/(isign*mmax)
wpr=-2.d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.d0
wi=0.d0
do m=1,mmax,2
do i=m,n,istep
j=i+mmax
tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
enddo
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
enddo
mmax=istep
goto 2
endif
return
END


After compiling:

[cva@ninja Oroginal]$ ifc tfourier.f
program TFOURIER
external subroutine REALFT
external subroutine FOUR1

123 Lines Compiled
[cva@ninja Oroginal]$ ./a.out
11.55388
-4.768372E-07
-2.384564
0.4743197
-0.2692006
0.1115067
-0.1038520
6.939162E-02
-4.609013E-02
4.609013E-02
-2.051385E-02
3.070122E-02
-7.870610E-03
1.900150E-02
-1.813536E-03
9.117555E-03