Beginner--How do I write output of 1-D Fourier Transform program? 1

[lilsalsa74]

I am a complete beginner at programming and using Fortran. I am actually a Physics major and I am trying to learn more about Fourier Transforms. I have been looking at Ooura's 1-D DFT program and it compiles but does not give an output of the Fourier Transform. How would I do this?
(Also, Dr. Ooura said on his site that it was OK to copy and use his program. Besides, I am merely using this as an educational tool.)
Here is the program:

! Fast Fourier/Cosine/Sine Transform
! dimension ne
! data length ower of 2
! decimation :frequency
! radix :2
! data :inplace
! table :not use
! -------- Complex DFT (Discrete Fourier Transform) --------
! [definition]
! <case1>
! X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n
! <case2>
! X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n
! (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
! [usage]
! <case1>
! call cdft(2*n, cos(pi/n), sin(pi/n), a)
! <case2>
! call cdft(2*n, cos(pi/n), -sin(pi/n), a)
! [parameters]
! 2*n :data length (integer)
! n >= 1, n = power of 2
! a(0:2*n-1) :input/output data (real*8)
! input data
! a(2*j) = Re(x(j)), a(2*j+1) = Im(x(j)), 0<=j<n
! output data
! a(2*k) = Re(X(k)), a(2*k+1) = Im(X(k)), 0<=k<n
! [remark]
! Inverse of
! call cdft(2*n, cos(pi/n), -sin(pi/n), a)
! is
! call cdft(2*n, cos(pi/n), sin(pi/n), a)
! do j = 0, 2 * n - 1
! a(j) = a(j) / n
! end do
subroutine cdft(n, wr, wi, a)
integer n, i, j, k, l, m
real*8 wr, wi, a(0 : n - 1), wmr, wmi, wkr, wki, wdr, wdi,
& ss, xr, xi
wmr = wr
wmi = wi
m = n
do while (m .gt. 4)
l = m / 2
wkr = 1
wki = 0
wdr = 1 - 2 * wmi * wmi
wdi = 2 * wmi * wmr
ss = 2 * wdi
wmr = wdr
wmi = wdi
do j = 0, n - m, m
i = j + l
xr = a(j) - a(i)
xi = a(j + 1) - a(i + 1)
a(j) = a(j) + a(i)
a(j + 1) = a(j + 1) + a(i + 1)
a(i) = xr
a(i + 1) = xi
xr = a(j + 2) - a(i + 2)
xi = a(j + 3) - a(i + 3)
a(j + 2) = a(j + 2) + a(i + 2)
a(j + 3) = a(j + 3) + a(i + 3)
a(i + 2) = wdr * xr - wdi * xi
a(i + 3) = wdr * xi + wdi * xr
end do
do k = 4, l - 4, 4
wkr = wkr - ss * wdi
wki = wki + ss * wdr
wdr = wdr - ss * wki
wdi = wdi + ss * wkr
do j = k, n - m + k, m
i = j + l
xr = a(j) - a(i)
xi = a(j + 1) - a(i + 1)
a(j) = a(j) + a(i)
a(j + 1) = a(j + 1) + a(i + 1)
a(i) = wkr * xr - wki * xi
a(i + 1) = wkr * xi + wki * xr
xr = a(j + 2) - a(i + 2)
xi = a(j + 3) - a(i + 3)
a(j + 2) = a(j + 2) + a(i + 2)
a(j + 3) = a(j + 3) + a(i + 3)
a(i + 2) = wdr * xr - wdi * xi
a(i + 3) = wdr * xi + wdi * xr
end do
end do
m = l
end do
if (m .gt. 2) then
do j = 0, n - 4, 4
xr = a(j) - a(j + 2)
xi = a(j + 1) - a(j + 3)
a(j) = a(j) + a(j + 2)
a(j + 1) = a(j + 1) + a(j + 3)
a(j + 2) = xr
a(j + 3) = xi
end do
end if
if (n .gt. 4) call bitrv2(n, a)
end
subroutine bitrv2(n, a)
integer n, j, j1, k, k1, l, m, m2, n2
real*8 a(0 : n - 1), xr, xi
m = n / 4
m2 = 2 * m
n2 = n - 2
k = 0
do j = 0, m2 - 4, 4
if (j .lt. k) then
xr = a(j)
xi = a(j + 1)
a(j) = a(k)
a(j + 1) = a(k + 1)
a(k) = xr
a(k + 1) = xi
else if (j .gt. k) then
j1 = n2 - j
k1 = n2 - k
xr = a(j1)
xi = a(j1 + 1)
a(j1) = a(k1)
a(j1 + 1) = a(k1 + 1)
a(k1) = xr
a(k1 + 1) = xi
end if
k1 = m2 + k
xr = a(j + 2)
xi = a(j + 3)
a(j + 2) = a(k1)
a(j + 3) = a(k1 + 1)
a(k1) = xr
a(k1 + 1) = xi
l = m
do while (k .ge. l)
k = k - l
l = l / 2
end do
k = k + l
end do
end

 

[mikrom]

What you posted are only subroutines. The subroutine bitrv2() is needed for the subroutine cdft().
You need to write the main program, which calls the subroutine cdft().

 

[lilsalsa74]

So if I wrote a main program:

program main
call cdft(2*n, cos(pi/n), sin(pi/n), a)

How would I go about inputting data and then outputting the result of the fourier transform? Would I create an array a?
Any help would be much appreciated!!

 

[lilsalsa74]

also, thank you for your reply mikrom!

 

[mikrom]

You should first create the array a, then call the procedure cdft() and then write the result. But I don't know what is the result of the subroutine - it seems that the procedure modifies the array a and that is the result.

 

[lilsalsa74]

Ok so I tried creating a main program with an array a(i). I then tried calling a different subroutine which computes a Real Discrete Fourier Transform since I gave a(i) real values. Here is the main program:

program maincall
real a(4)
real, parameter :: pi = 3.14159265
integer n
n = 4
a(1) = -1.0
a(2) = 1.0
a(3) = -1.0
a(4) = 1.0
write(*,*)'The original data:',a
call rdft(n, cos(pi/n), sin(pi/n), a)
write(*,*)'The Fourier coefficients',a
end

Here is the RDFT subroutine:
! -------- Real DFT / Inverse of Real DFT --------
! [definition]
! <case1> RDFT
! R(k) = sum_j=0^n-1 a(j)*cos(2*pi*j*k/n), 0<=k<=n/2
! I(k) = sum_j=0^n-1 a(j)*sin(2*pi*j*k/n), 0<k<n/2
! <case2> IRDFT (excluding scale)
! a(k) = R(0)/2 + R(n/2)/2 +
! sum_j=1^n/2-1 R(j)*cos(2*pi*j*k/n) +
! sum_j=1^n/2-1 I(j)*sin(2*pi*j*k/n), 0<=k<n
! [usage]
! <case1>
! call rdft(n, cos(pi/n), sin(pi/n), a)
! <case2>
! call rdft(n, cos(pi/n), -sin(pi/n), a)
! [parameters]
! n :data length (integer)
! n >= 2, n = power of 2
! a(0:n-1) :input/output data (real*8)
! <case1>
! output data
! a(2*k) = R(k), 0<=k<n/2
! a(2*k+1) = I(k), 0<k<n/2
! a(1) = R(n/2)
! <case2>
! input data
! a(2*j) = R(j), 0<=j<n/2
! a(2*j+1) = I(j), 0<j<n/2
! a(1) = R(n/2)
! [remark]
! Inverse of
! call rdft(n, cos(pi/n), sin(pi/n), a)
! is
! call rdft(n, cos(pi/n), -sin(pi/n), a)
! do j = 0, n - 1
! a(j) = a(j) * 2 / n
! end do
subroutine rdft(n, wr, wi, a)
integer n, j, k
real*8 wr, wi, a(0 : n - 1), wmr, wmi, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi
if (n .gt. 4) then
wkr = 0
wki = 0
wdr = wi * wi
wdi = wi * wr
ss = 4 * wdi
wmr = 1 - 2 * wdr
wmi = 2 * wdi
if (wmi .ge. 0) then
call cdft(n, wmr, wmi, a)
xi = a(0) - a(1)
a(0) = a(0) + a(1)
a(1) = xi
end if
do k = n / 2 - 4, 4, -4
j = n - k
xr = a(k + 2) - a(j - 2)
xi = a(k + 3) + a(j - 1)
yr = wdr * xr - wdi * xi
yi = wdr * xi + wdi * xr
a(k + 2) = a(k + 2) - yr
a(k + 3) = a(k + 3) - yi
a(j - 2) = a(j - 2) + yr
a(j - 1) = a(j - 1) - yi
wkr = wkr + ss * wdi
wki = wki + ss * (0.5d0 - wdr)
xr = a(k) - a(j)
xi = a(k + 1) + a(j + 1)
yr = wkr * xr - wki * xi
yi = wkr * xi + wki * xr
a(k) = a(k) - yr
a(k + 1) = a(k + 1) - yi
a(j) = a(j) + yr
a(j + 1) = a(j + 1) - yi
wdr = wdr + ss * wki
wdi = wdi + ss * (0.5d0 - wkr)
end do
j = n - 2
xr = a(2) - a(j)
xi = a(3) + a(j + 1)
yr = wdr * xr - wdi * xi
yi = wdr * xi + wdi * xr
a(2) = a(2) - yr
a(3) = a(3) - yi
a(j) = a(j) + yr
a(j + 1) = a(j + 1) - yi
if (wmi .lt. 0) then
a(1) = 0.5d0 * (a(0) - a(1))
a(0) = a(0) - a(1)
call cdft(n, wmr, wmi, a)
end if
else
if (wi .lt. 0) then
a(1) = 0.5d0 * (a(0) - a(1))
a(0) = a(0) - a(1)
end if
if (n .gt. 2) then
xr = a(0) - a(2)
xi = a(1) - a(3)
a(0) = a(0) + a(2)
a(1) = a(1) + a(3)
a(2) = xr
a(3) = xi
end if
if (wi .ge. 0) then
xi = a(0) - a(1)
a(0) = a(0) + a(1)
a(1) = xi
end if
end if
end

I tried to run the program and I got the following run-time error:

Attempt to call a routine with argument number two as a real(kind=1) when a real(kind=2) was required.

I'm a bit confused. I did, at least, get the array to print out before it had been called by the RDFT program.

 

[lilsalsa74]

Thank you so much again, Mikrom. Just having someone give input really helps a lot.
So I changed the declaration of pi and a(4) to kind=2 and now the program does, in fact, run. The values of the coefficients the program outputs seem to be correct(I checked with a well-known function) but they are in the wrong order. I will look more into this.
If you have any more suggestions, please feel free to let me know. Thank you again!!!

 

[mikrom]

I named the file with your subroutine to
fft.f90

[COLOR=#0000ff]! Fast Fourier/Cosine/Sine Transform[/color]
[COLOR=#0000ff]!     dimension   :one[/color]
[COLOR=#0000ff]!     data length :power of 2[/color]
[COLOR=#0000ff]!     decimation  :frequency[/color]
[COLOR=#0000ff]!     radix       :2[/color]
[COLOR=#0000ff]!     data        :inplace[/color]
[COLOR=#0000ff]!     table       :not use[/color]
[COLOR=#0000ff]! -------- Complex DFT (Discrete Fourier Transform) --------[/color]
[COLOR=#0000ff]!     [definition][/color]
[COLOR=#0000ff]!         <case1>[/color]
[COLOR=#0000ff]!             X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n[/color]
[COLOR=#0000ff]!         <case2>[/color]
[COLOR=#0000ff]!             X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n[/color]
[COLOR=#0000ff]!         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)[/color]
[COLOR=#0000ff]!     [usage][/color]
[COLOR=#0000ff]!         <case1>[/color]
[COLOR=#0000ff]!             call cdft(2*n, cos(pi/n), sin(pi/n), a)[/color]
[COLOR=#0000ff]!         <case2>[/color]
[COLOR=#0000ff]!             call cdft(2*n, cos(pi/n), -sin(pi/n), a)[/color]
[COLOR=#0000ff]!     [parameters][/color]
[COLOR=#0000ff]!         2*n        :data length (integer)[/color]
[COLOR=#0000ff]!                     n >= 1, n = power of 2[/color]
[COLOR=#0000ff]!         a(0:2*n-1) :input/output data (real*8)[/color]
[COLOR=#0000ff]!                     input data[/color]
[COLOR=#0000ff]!                         a(2*j) = Re(x(j)), a(2*j+1) = Im(x(j)), 0<=j<n[/color]
[COLOR=#0000ff]!                     output data[/color]
[COLOR=#0000ff]!                         a(2*k) = Re(X(k)), a(2*k+1) = Im(X(k)), 0<=k<n[/color]
[COLOR=#0000ff]!     [remark][/color]
[COLOR=#0000ff]!         Inverse of [/color]
[COLOR=#0000ff]!             call cdft(2*n, cos(pi/n), -sin(pi/n), a)[/color]
[COLOR=#0000ff]!         is [/color]
[COLOR=#0000ff]!             call cdft(2*n, cos(pi/n), sin(pi/n), a)[/color]
[COLOR=#0000ff]!             do j = 0, 2 * n - 1[/color]
[COLOR=#0000ff]!                 a(j) = a(j) / n[/color]
[COLOR=#0000ff]!             end do[/color]
[COLOR=#a020f0]subroutine[/color] cdft(n, wr, wi, a)
      [COLOR=#2e8b57][b]integer[/b][/color] n, i, j, k, l, m
[COLOR=#2e8b57][b]      real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] wr, wi, a([COLOR=#ff00ff]0[/color] : n [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]1[/color]), wmr, wmi, wkr, wki, wdr, wdi, [highlight #ffff00][COLOR=#0000ff]&[/color][/highlight]
             ss, xr, xi
      wmr [COLOR=#804040][b]=[/b][/color] wr
      wmi [COLOR=#804040][b]=[/b][/color] wi
      m [COLOR=#804040][b]=[/b][/color] n
      [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] (m [COLOR=#804040][b].gt.[/b][/color] [COLOR=#ff00ff]4[/color])
          l [COLOR=#804040][b]=[/b][/color] m [COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2[/color]
          wkr [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
          wki [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
          wdr [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color] [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#804040][b]*[/b][/color] wmi [COLOR=#804040][b]*[/b][/color] wmi
          wdi [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#804040][b]*[/b][/color] wmi [COLOR=#804040][b]*[/b][/color] wmr
          ss [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#804040][b]*[/b][/color] wdi
          wmr [COLOR=#804040][b]=[/b][/color] wdr
          wmi [COLOR=#804040][b]=[/b][/color] wdi
          [COLOR=#804040][b]do[/b][/color] j [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color], n [COLOR=#804040][b]-[/b][/color] m, m
              i [COLOR=#804040][b]=[/b][/color] j [COLOR=#804040][b]+[/b][/color] l
              xr [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]-[/b][/color] a(i)
              xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(j) [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]+[/b][/color] a(i)
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(i) [COLOR=#804040][b]=[/b][/color] xr
              a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] xi
              xr [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
              xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
              a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] wdr [COLOR=#804040][b]*[/b][/color] xr [COLOR=#804040][b]-[/b][/color] wdi [COLOR=#804040][b]*[/b][/color] xi
              a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] wdr [COLOR=#804040][b]*[/b][/color] xi [COLOR=#804040][b]+[/b][/color] wdi [COLOR=#804040][b]*[/b][/color] xr
          [COLOR=#804040][b]end do[/b][/color]
          [COLOR=#804040][b]do[/b][/color] k [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]4[/color], l [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]4[/color]
              wkr [COLOR=#804040][b]=[/b][/color] wkr [COLOR=#804040][b]-[/b][/color] ss [COLOR=#804040][b]*[/b][/color] wdi
              wki [COLOR=#804040][b]=[/b][/color] wki [COLOR=#804040][b]+[/b][/color] ss [COLOR=#804040][b]*[/b][/color] wdr
              wdr [COLOR=#804040][b]=[/b][/color] wdr [COLOR=#804040][b]-[/b][/color] ss [COLOR=#804040][b]*[/b][/color] wki
              wdi [COLOR=#804040][b]=[/b][/color] wdi [COLOR=#804040][b]+[/b][/color] ss [COLOR=#804040][b]*[/b][/color] wkr
              [COLOR=#804040][b]do[/b][/color] j [COLOR=#804040][b]=[/b][/color] k, n [COLOR=#804040][b]-[/b][/color] m [COLOR=#804040][b]+[/b][/color] k, m
                  i [COLOR=#804040][b]=[/b][/color] j [COLOR=#804040][b]+[/b][/color] l
                  xr [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]-[/b][/color] a(i)
                  xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
                  a(j) [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]+[/b][/color] a(i)
                  a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
                  a(i) [COLOR=#804040][b]=[/b][/color] wkr [COLOR=#804040][b]*[/b][/color] xr [COLOR=#804040][b]-[/b][/color] wki [COLOR=#804040][b]*[/b][/color] xi
                  a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] wkr [COLOR=#804040][b]*[/b][/color] xi [COLOR=#804040][b]+[/b][/color] wki [COLOR=#804040][b]*[/b][/color] xr
                  xr [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
                  xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]-[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
                  a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
                  a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]+[/b][/color] a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
                  a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] wdr [COLOR=#804040][b]*[/b][/color] xr [COLOR=#804040][b]-[/b][/color] wdi [COLOR=#804040][b]*[/b][/color] xi
                  a(i [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] wdr [COLOR=#804040][b]*[/b][/color] xi [COLOR=#804040][b]+[/b][/color] wdi [COLOR=#804040][b]*[/b][/color] xr
              [COLOR=#804040][b]end do[/b][/color]
          [COLOR=#804040][b]end do[/b][/color]
          m [COLOR=#804040][b]=[/b][/color] l
      [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#804040][b]if[/b][/color] (m [COLOR=#804040][b].gt.[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]then[/b][/color]
          [COLOR=#804040][b]do[/b][/color] j [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color], n [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]4[/color]
              xr [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]-[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
              xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]-[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
              a(j) [COLOR=#804040][b]=[/b][/color] a(j) [COLOR=#804040][b]+[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] xr
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] xi
          [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#804040][b]end if[/b][/color]
      [COLOR=#804040][b]if[/b][/color] (n [COLOR=#804040][b].gt.[/b][/color] [COLOR=#ff00ff]4[/color]) [COLOR=#a020f0]call[/color] bitrv2(n, a)
      [COLOR=#a020f0]end[/color]
[COLOR=#a020f0]subroutine[/color] bitrv2(n, a)
      [COLOR=#2e8b57][b]integer[/b][/color] n, j, j1, k, k1, l, m, m2, n2
[COLOR=#2e8b57][b]      real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] a([COLOR=#ff00ff]0[/color] : n [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]1[/color]), xr, xi
      m [COLOR=#804040][b]=[/b][/color] n [COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]4[/color]
      m2 [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#804040][b]*[/b][/color] m
      n2 [COLOR=#804040][b]=[/b][/color] n [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]2[/color]
      k [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
      [COLOR=#804040][b]do[/b][/color] j [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color], m2 [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]4[/color]
          [COLOR=#804040][b]if[/b][/color] (j [COLOR=#804040][b].lt.[/b][/color] k) [COLOR=#804040][b]then[/b][/color]
              xr [COLOR=#804040][b]=[/b][/color] a(j)
              xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(j) [COLOR=#804040][b]=[/b][/color] a(k)
              a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] a(k [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(k) [COLOR=#804040][b]=[/b][/color] xr
              a(k [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] xi
          [COLOR=#804040][b]else[/b][/color] [COLOR=#804040][b]if[/b][/color] (j [COLOR=#804040][b].gt.[/b][/color] k) [COLOR=#804040][b]then[/b][/color]
              j1 [COLOR=#804040][b]=[/b][/color] n2 [COLOR=#804040][b]-[/b][/color] j
              k1 [COLOR=#804040][b]=[/b][/color] n2 [COLOR=#804040][b]-[/b][/color] k
              xr [COLOR=#804040][b]=[/b][/color] a(j1)
              xi [COLOR=#804040][b]=[/b][/color] a(j1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(j1) [COLOR=#804040][b]=[/b][/color] a(k1)
              a(j1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] a(k1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
              a(k1) [COLOR=#804040][b]=[/b][/color] xr
              a(k1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] xi
          [COLOR=#804040][b]end if[/b][/color]
          k1 [COLOR=#804040][b]=[/b][/color] m2 [COLOR=#804040][b]+[/b][/color] k
          xr [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color])
          xi [COLOR=#804040][b]=[/b][/color] a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color])
          a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] a(k1)
          a(j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] a(k1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color])
          a(k1) [COLOR=#804040][b]=[/b][/color] xr
          a(k1 [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] xi
          l [COLOR=#804040][b]=[/b][/color] m
          [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] (k [COLOR=#804040][b].ge.[/b][/color] l)
              k [COLOR=#804040][b]=[/b][/color] k [COLOR=#804040][b]-[/b][/color] l
              l [COLOR=#804040][b]=[/b][/color] l [COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2[/color]
          [COLOR=#804040][b]end do[/b][/color]
          k [COLOR=#804040][b]=[/b][/color] k [COLOR=#804040][b]+[/b][/color] l
      [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#a020f0]end[/color]

and then I created the main program
fft_example.f90

[COLOR=#a020f0]program[/color] fft_example
  [COLOR=#2e8b57][b]integer[/b][/color] :: n
[COLOR=#2e8b57][b]  real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] :: a([COLOR=#ff00ff]0[/color]:[COLOR=#ff00ff]3[/color])

  [COLOR=#0000ff]! enter input values[/color]
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Input values:'[/color]
  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]4[/color]
  a([COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color]
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'a['[/color],i,[COLOR=#ff00ff]']= '[/color], a(i)
  [COLOR=#804040][b]end do[/b][/color]

  [COLOR=#0000ff]! call FFT procedure[/color]
  [COLOR=#a020f0]call[/color] cdft(n, [COLOR=#008080]dcos[/color]([COLOR=#008080]dfloat[/color](pi[COLOR=#804040][b]/[/b][/color]n)), [COLOR=#008080]dsin[/color]([COLOR=#008080]dfloat[/color](pi[COLOR=#804040][b]/[/b][/color]n)), a)

  [COLOR=#0000ff]! print output values[/color]
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Output values:'[/color]  
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'a['[/color],i,[COLOR=#ff00ff]']= '[/color], a(i)
  [COLOR=#804040][b]end do[/b][/color]
[COLOR=#a020f0]end program[/color] fft_example

[COLOR=#0000ff]! include source with FFT subroutine[/color]
[COLOR=#a020f0]include[/color] [COLOR=#ff00ff]'fft.f90'[/color]

Compilation and running

$ g95 fft_example.f90 -o fft_example

$ fft_example
 Input values:
 a[ 0 ]=  -1.
 a[ 1 ]=  1.
 a[ 2 ]=  -1.
 a[ 3 ]=  1.
 Output values:
 a[ 0 ]=  -2.
 a[ 1 ]=  2.
 a[ 2 ]=  0.
 a[ 3 ]=  0.

 

[mikrom]

I have an error in my program. I thought when Fortran is for mathematics it knows the value of PI implicitely
But it's not true - the value of PI must be given (or computed) explicitely.

So here is a the corected source
fft_example.f90

[COLOR=#a020f0]program[/color] fft_example
  [COLOR=#2e8b57][b]integer[/b][/color] :: n
[COLOR=#2e8b57][b]  real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] :: a([COLOR=#ff00ff]0[/color]:[COLOR=#ff00ff]3[/color]), pi

  [COLOR=#0000ff]! enter input values[/color]
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Input values:'[/color]
  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]4[/color]
  a([COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.0[/color]
  a([COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color]
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'a['[/color],i,[COLOR=#ff00ff]']= '[/color], a(i)
  [COLOR=#804040][b]end do[/b][/color]

  [COLOR=#0000ff]! compute pi[/color]
  pi [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color] [COLOR=#804040][b]*[/b][/color] [COLOR=#008080]DACOS[/color]([COLOR=#008080]dfloat[/color]([COLOR=#ff00ff]0.0[/color]))
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Pi = '[/color], pi
  [COLOR=#0000ff]! call FFT procedure[/color]
  [COLOR=#a020f0]call[/color] cdft(n, [COLOR=#008080]dcos[/color](pi[COLOR=#804040][b]/[/b][/color]n), [COLOR=#008080]dsin[/color](pi[COLOR=#804040][b]/[/b][/color]n), a)

  [COLOR=#0000ff]! print output values[/color]
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Output values:'[/color]  
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'a['[/color],i,[COLOR=#ff00ff]']= '[/color], a(i)
  [COLOR=#804040][b]end do[/b][/color]
[COLOR=#a020f0]end program[/color] fft_example

[COLOR=#0000ff]! include source with FFT subroutine[/color]
[COLOR=#a020f0]include[/color] [COLOR=#ff00ff]'fft.f90'[/color]

The intrinsic function dfloat() converts single precision float (i.e. REAL*4) to double precision (i.e. REAL*8)

Output:

$ g95 fft_example.f90 -o fft_example

$ fft_example
 Input values:
 a[ 0 ]=  -1.
 a[ 1 ]=  1.
 a[ 2 ]=  -1.
 a[ 3 ]=  1.
 Pi =  3.141592653589793
 Output values:
 a[ 0 ]=  -2.
 a[ 1 ]=  2.
 a[ 2 ]=  0.
 a[ 3 ]=  0.

 

[lilsalsa74]

Thank you so much!! That was extremely helpful. Aren't Fourier Transforms fun?!

 

[mikrom]

YES !