If all you need is a simple

Pramudito__Ary · ‎10-29-2017

Hello everyone I am fortran enthusiast since compaq, but no so well. I have a got problem in fft. I see a lot of fft subroutine, and I don't know which is the correct one for my project. I got ambient acceleration time series real*8 (m/s^2 versus seconds) data from accelerometer of the structure. And I want to find the fundamental or eigen frequency using fft power spectra. As I broWsing, maybe it's look like matlab can done in https://www.mathworks.com/help/matlab/examples/fft-for-spectral-analysis.html . So it's need fft and multiply by its conjugate. What function and subroutine that fit with my data. There is cosine , 1d/2d fft subroutine in web that I don't know what is for.

JVanB · ‎10-30-2017

Looks like you need 1D real to half complex FFT. You know, MKL comes with libraries to do this stuff, so you might try inquiring on the MKL forum.

dboggs · ‎11-01-2017

If all you need is a simple rudimentary fft, then getting MKL to work with your code may be way overkill unless you already have experience with it. An fft routine for this is very short; you can easily type it into your code way faster than you can get MKL to work. Here's the routine that I use, it is based on the code published by Newland in his book on random vibration.

(Sorry about the formatting here; I just can't get a handle on it)

Once you have the complex fft results in {Z}, then get the PSD by using the complex multiplication

<pre class="brush:Fortran">
      DO iz = 1, nft/2 ! OK to store psd values in first half of {z}
         zreal = real (z(iz))
         zimag = aimag (z(iz))
         z(iz) = cmplx (2. * (zreal**2 + zimag**2), 0.)
      END DO
</pre>

<pre class="brush:Fortran">!===============================================================================
! FFT1 (Newland)
!-------------------------------------------------------------------------------
! FFT routine given by Newland, p.220
!
! Called by: PSD_CALC
!
! Calls:      None
!
! In:
!     {Z} -    Vector of complex points xi, i = 1,...NB
!     N -      Log2 of the number of points
!     NB -     No. of points, = 2**N
!
! Out:
!     {Z} -    Vector of complex Fourier transform Xi, i = 1, 2, ...NB.
!              i = 1 corresponds to frequency 0
!              i = NB/2+1 corresponds to Nyquist frequency
!              i > NB/2+1 are symmetrically redundant.
!     Example: NB = 16, N = 4
!     X1 corresponds to frequency 0
!     X2 corresponds to frequency 1
!     X8 corresponds to frequency 7
!     X9 corresponds t

!===============================================================================
!  FFT1 (Newland)
!-------------------------------------------------------------------------------
!  FFT routine given by Newland, p.220
!
!  Called by:  PSD_CALC
!
!  Calls:      None
!
!  In:
!     {Z} -    Vector of complex points xi, i = 1,...NB
!     N -      Log2 of the number of points
!     NB -     No. of points, = 2**N
!
!  Out:
!     {Z} -    Vector of complex Fourier transform Xi, i = 1, 2, ...NB.
!              i = 1 corresponds to frequency 0
!              i = NB/2+1 corresponds to Nyquist frequency
!              i > NB/2+1 are symmetrically redundant.
!     Example: NB = 16, N = 4
!     X1  corresponds to frequency 0
!     X2  corresponds to frequency 1
!     X8  corresponds to frequency 7
!     X9  corresponds to frequency 8 = Nyquist frequency
!     X10 corresponds to frequency 9 = X8
!     X11 corresponds to frequency 10 = X7
!     X16 corresponds to frequency 15 = X1
!
!     These results are more sensical if Z is indexed to begin with 0 instead of
!     1. It is left at 1 here for conformance with Newland. However, the calling
!     routine may use Z(0:NB-1). 
!-------------------------------------------------------------------------------
   SUBROUTINE fft_Newland (Z, N, NB)
   
      IMPLICIT NONE
     ! Arguments
      INTEGER(4)  :: N, NB
      COMPLEX(8)     :: Z(NB)
     ! Locals
      INTEGER(4)  :: NBD2, NBM1
      COMPLEX(8)     :: U, W, T
      INTEGER(4)  :: J, K, L, M, ME, LPK
      REAL(8), PARAMETER   :: PI = 3.14159265
      real(8)     :: rnb
      
      rnb = nb
      DO  J = 1, NB
          Z(J) = Z(J) / rNB
      END DO

!  Reorder sequence according to Fig. 12.8.

      NBD2 = NB / 2
      NBM1 = NB - 1
      J = 1
      DO L = 1, NBM1
         IF (l < J) THEN
            T = Z(J)
            Z(J) = Z(L)
            Z(L) = T
         END IF
         K = NBD2
 30      IF (K >= J) GO TO 40
         J = J - K
         K = K / 2
         GO TO 30
 40      J = J + K
      END DO ! next L

!  Calculate FFT according to Fig. 12.5

      DO M = 1, N
         U = (1.0, 0.0)
         ME = 2**M
         K = ME / 2
         W = CMPLX (COS(PI/K), -SIN(PI/K))
         DO J = 1, K
            DO L = J, NB, ME
               LPK = L + K
               T = Z(LPK) * U
               Z(LPK) = Z(L) - T
               Z(L) = Z(L) + T
            END DO ! L
            U = U * W
         END DO ! J
      END DO ! M

      RETURN
   END SUBROUTINE fft_Newland</pre>

Pramudito__Ary · ‎11-08-2017

Thank you for respond

Ok i had the response for accelerometer with 1921 data for 60 seconds. So my X array is about X(1921) and Y(1921). X is the time second data, Y is the amplitudo.

To get fft one array real that i need a 1D real to complex FFT. I saw maybe fftw3 is the answer. (any mkl fo, or even cxml / imsl compaq info , i will like to hear it)

I compile the 32 bit fortran wrapper library via administrator command prompt at mkl\interfaces\fftw2xf using nmake lib32

I put the library, fft3.fi and fftw.h from intel mkl\include\fftw to the workspace.

Somehow it need declaration integer FFTW_ESTIMATE, i dont know where to find FFTW_ESTIMATE at already header files, so i put it my subroutine.

Question --

From complex to find the Power spectrum to find it...

1. How connect to call newland FFT dboogs, seems it need an 'even data'.

2. I saw FFTw2 examples in C, how implemented in Fortran using fftw3 ? Here also we need only half+1 data complex.

3. What i need for complex data var: full data num + 1, or half + 1 ?, how to correspond it to Hz ?

4. How to multiply with conjugate like in above matlab command : > real Pyy = Y.*conj(Y)/totalsampling; where the Y is the complex data FFT result. What conjugate function in fortran?

Regards

subroutine fftwhat

     implicit none
 

!     The fftw variable that needed to paste...where the header files?? 
      integer FFTW_FORWARD,FFTW_BACKWARD
      parameter (FFTW_FORWARD=-1,FFTW_BACKWARD=1)

      integer FFTW_REAL_TO_COMPLEX,FFTW_COMPLEX_TO_REAL
      parameter (FFTW_REAL_TO_COMPLEX=-1,FFTW_COMPLEX_TO_REAL=1)

      integer FFTW_ESTIMATE,FFTW_MEASURE
      parameter (FFTW_ESTIMATE=0,FFTW_MEASURE=1)

      integer FFTW_OUT_OF_PLACE,FFTW_IN_PLACE,FFTW_USE_WISDOM
      parameter (FFTW_OUT_OF_PLACE=0)
      parameter (FFTW_IN_PLACE=8,FFTW_USE_WISDOM=16)

      integer FFTW_THREADSAFE
      parameter (FFTW_THREADSAFE=128)

!     End of variable fftw3 header

     REAL*8, ALLOCATABLE :: x(:)
     REAL*8, ALLOCATABLE :: y(:)
     REAL*8, ALLOCATABLE :: y2(:)
     complex*8, ALLOCATABLE :: coef(:)
     integer totalsampling,i

     integer  plan_forward

 
     totalsampling = 1921
     
 
     allocate (x(totalsampling)) !time
     allocate (y(totalsampling)) !amplitudo
     allocate (coef(totalsampling+1)) !for fft   DO i need half+1 or total sampling + 1 ??

 
      do i=1,totalsampling

        x(i) = 0.d0 + (0.31250E-01* (i-1))

     enddo

    do i=1,totalsampling

        y(i) = sin(i*3.14/10)

     enddo

     call dfftw_plan_dft_r2c_1d ( plan_forward, totalsampling, y, coef, FFTW_ESTIMATE )

    call dfftw_execute ( plan_forward )

!  From here i got the complex coef right? 

! How to connect with above dboogs newland fft? or power spectra



!  Release the memory associated with the plans.
     call dfftw_destroy_plan ( plan_forward )


end

The C examples spectra...

! http://csweb.cs.wfu.edu/~torgerse/fftw_2.1.2/fftw_2.html

#include <rfftw.h>
...
{
     fftw_real in, out, power_spectrum[N/2+1];
     rfftw_plan p;
     int k;
     ...
     p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
     ...
     rfftw_one(p, in, out);
     power_spectrum[0] = out[0]*out[0];  /* DC component */
     for (k = 1; k < (N+1)/2; ++k)  /* (k < N/2 rounded up) */
          power_spectrum = out*out + out[N-k]*out[N-k];
     if (N % 2 == 0) /* N is even */
          power_spectrum[N/2] = out[N/2]*out[N/2];  /* Nyquist freq. */
     ...
     rfftw_destroy_plan(p);
}

JVanB · ‎11-08-2017

N = 1921 = 17*113 = (2**4+1)*(7*2**4+1) is fairly efficient via the PFA+Rader methods but I haven't seen an efficient Rader FFT of order 113 and it takes a day or more to write those things out by hand. The alternative is the chirp-Z fft and this is how most packages would probably handle it because all it needs to work is a power of 2 fft. I have seen the claim that it's always straightforward to make a corresponding real-half complex fft from any full complex fft but I have never seen the claimants back up their statements in the case of a chirp-Z fft! Even so, the factor of 8 for using full complex chirp-Z instead of real-half complex PFA+Rader will probably be negligible for such small data sets. I would try to make this fly with MKL rather that FFTW if at all possible.

Pramudito__Ary · ‎11-10-2017

Thank you for dboogs and Repeatoffender.

After reading vibrationdata.wordpress.com , i realized 2^n data is more easier to transform. So data will be truncated to 2^n sampling. There is also PSD and FFT Fortran there that i check today.

For dboogs, the newland FFT is stuck on ordering when i try to verify yesterday with sin signal.

program example_fftw
implicit none

include 'C:\Program Files (x86)\Intel\Composer XE 2015\mkl\include\fftw\fftw3.f'
!*****************************************************************************80
!
!! TEST02: real 1D data.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2011
!
!  Author:
!
!    John Burkardt
!

  

  integer ( kind = 4 ), parameter :: n = 256
  integer ( kind = 4 ), parameter :: nc = 129  
  character (len = 256 ) csv_file_name
  integer ( kind = 4 ) i
  real ( kind = 8 ) in(n)
  real ( kind = 8 ) in2(n)
  complex ( kind = 8 ) out(nc)
  integer ( kind = 8 ) plan_backward
  integer ( kind = 8 ) plan_forward
  real*8 samplesec, sec(256),pi,secakhir
  
  csv_file_name = "datainput_fftw.log"
  open(unit=20,file=csv_file_name,status='replace',access ='sequential',form='formatted',action='write')  

  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) 'TEST02'
  write ( 20, '(a)' ) '  Demonstrate FFTW3 on a single vector'
  write ( 20, '(a)' ) '  of real data.'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Transform data to FFT coefficients.'
  write ( 20, '(a)' ) '  Backtransform FFT coefficients to recover '
  write ( 20, '(a)' ) '  data.'
  write ( 20, '(a)' ) '  Compare recovered data to original data.'
!
!  Set up the input data, a real vector of length N.
!
!  call r8vec_uniform_01 ( n, seed, in )
!
  samplesec = 0.001d0
  pi = acos(dble(-1.0))
  do i=1,n
      sec(i) = 0.d0 + samplesec*(i-1)
      secakhir = sec(i)
      !in(i) = sin(2*pi*50*sec(i)) + sin(2*pi*120*sec(i));
      in(i) = sin(2*pi*50*sec(i))
  enddo
 
  
  
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Input Data:'
  write ( 20, '(a)' ) ' '

  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, in(i)
  end do
!
!  Set up a plan, and execute the plan to transform the IN data to
!  the OUT FFT coefficients.
!
  
  call dfftw_plan_dft_r2c_1d ( plan_forward, n, in, out, FFTW_ESTIMATE )

  call dfftw_execute ( plan_forward )
  write (*,*) 'FFT Forward done'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Output FFT Coefficients:'
  write ( 20, '(a)' ) ' '

  do i = 1, nc
    write ( 20, '(2x,i4,2x,2g14.6)' ) i, out(i)
  end do
!
!  Set up a plan, and execute the plan to backtransform the
!  complex FFT coefficients in OUT to real data.
!
  call dfftw_plan_dft_c2r_1d ( plan_backward, n, out, in2, FFTW_ESTIMATE )

  call dfftw_execute ( plan_backward )
  write (*,*) 'FFT Backward done'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Recovered input data divide by N:'
  write ( 20, '(a)' ) ' '

  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, in2(i) / real ( n, kind = 8 )
  end do
!
!  Release the memory associated with the plans.
!
  call dfftw_destroy_plan ( plan_forward )
  call dfftw_destroy_plan ( plan_backward )
  
  write ( 20,*) 'difference'
  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, (in2(i) / real ( n, kind = 8 )) - in(i) 
  end do
! Trying PSD
  call hey(out,n,nc)

  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Output Newland'
  write ( 20, '(a)' ) ' '

 do i = 1, nc
   write ( 20, '(2x,i4,2x,2g14.6)' ) i, out(i)
 end do

  
  close(20)
end 


subroutine hey(out,n,nc)
implicit none
  !integer ( kind = 4 ), parameter :: n = 256
  !integer ( kind = 4 ), parameter :: nc = 129  
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nc
  complex ( kind = 8 ) out(nc),out2(nc)
  real*8 llooog,llegup,llegdown,zreal,zimag
  integer iz,lleg
  
   llooog = dLOG (2.d0)
   llooog = log (real(n)) / llooog
   llegup = ceiling(llooog)
   llegdown = floor(llooog)
   if (abs(llegup - llooog).lt.abs(llegdown - llooog)) lleg = llegup
   if (abs(llegdown - llooog).lt.abs(llegup  - llooog)) lleg = llegdown
    
      DO iz = 1, nc ! OK to store psd values in first half of {z}
         zreal = real (out(iz))
         zimag = aimag (out(iz))
         out2(iz) = cmplx (2. * (zreal**2 + zimag**2), 0.)
      END DO
    
    
    
    !pause
     call fft_Newland(out2,  lleg, nc)

    end

       !===============================================================================

!  FFT1 (Newland)

!-------------------------------------------------------------------------------

!  FFT routine given by Newland, p.220

!

!  Called by:  PSD_CALC

!

!  Calls:      None

!

!  In:

!     {Z} -    Vector of complex points xi, i = 1,...NB

!     N -      Log2 of the number of points

!     NB -     No. of points, = 2**N

!

!  Out:

!     {Z} -    Vector of complex Fourier transform Xi, i = 1, 2, ...NB.

!              i = 1 corresponds to frequency 0

!              i = NB/2+1 corresponds to Nyquist frequency

!              i > NB/2+1 are symmetrically redundant.

!     Example: NB = 16, N = 4

!     X1  corresponds to frequency 0

!     X2  corresponds to frequency 1

!     X8  corresponds to frequency 7

!     X9  corresponds to frequency 8 = Nyquist frequency

!     X10 corresponds to frequency 9 = X8

!     X11 corresponds to frequency 10 = X7

!     X16 corresponds to frequency 15 = X1

!

!     These results are more sensical if Z is indexed to begin with 0 instead of

!     1. It is left at 1 here for conformance with Newland. However, the calling

!     routine may use Z(0:NB-1).
 
!-------------------------------------------------------------------------------


    
   SUBROUTINE fft_Newland (Z, N, NB)
!  FFT routine given by Newland, p.220
!
!  Called by:  PSD_CALC
!
!  Calls:      None
!
!  In:
!     {Z} -    Vector of complex points xi, i = 1,...NB
!     N -      Log2 of the number of points
!     NB -     No. of points, = 2**N
    

      IMPLICIT NONE

     ! Arguments

      INTEGER(4)  :: N, NB

      COMPLEX(8)     :: Z(NB)

     ! Locals

      INTEGER(4)  :: NBD2, NBM1

      COMPLEX(8)     :: U, W, T

      INTEGER(4)  :: J, K, L, M, ME, LPK

      REAL(8), PARAMETER   :: PI = 3.14159265

      real(8)     :: rnb

       

      rnb = nb

      DO  J = 1, NB

          Z(J) = Z(J) / rNB

      END DO

 

!  Reorder sequence according to Fig. 12.8.

 

      NBD2 = NB / 2

      NBM1 = NB - 1

      J = 1

      DO L = 1, NBM1
          
          write (*,*) 'Reorder Newland :',l,' / ', nbm1

         IF (l < J) THEN

            T = Z(J)

            Z(J) = Z(L)

            Z(L) = T

         END IF

         K = NBD2

30      IF (K >= J) GO TO 40

         J = J - K

         K = K / 2

         GO TO 30

 40      J = J + K

      END DO ! next L

 

!  Calculate FFT according to Fig. 12.5

 
      write (*,*) 'FFT newland starts'
      DO M = 1, N

          write (*,*) 'FFT Newland ',m,' / ', n
          U = (1.0, 0.0)

         ME = 2**M

         K = ME / 2

         W = CMPLX (COS(PI/K), -SIN(PI/K))

         DO J = 1, K

          write (*,*) 'FFT Newland sub ',j,' / ', k
          
             DO L = J, NB, ME

               LPK = L + K

               T = Z(LPK) * U

               Z(LPK) = Z(L) - T

               Z(L) = Z(L) + T

            END DO ! L

            U = U * W

         END DO ! J

      END DO ! M

 

      RETURN

    END SUBROUTINE !fft_Newland
    
    !--------------------------------------------------------------------------------

ITs stuck in ordering....

While trying MKL, i realize the output is complex on real variable. How to use it on complex input like on Newland FFT ?

program intelMKL1D


  use MKL_DFTI, forget => DFTI_DOUBLE, DFTI_DOUBLE => DFTI_DOUBLE_R
  implicit none
  ! Size of 1D transform
 
!*****************************************************************************80
!
!! TEST02: real 1D data.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2011
!
!  Author:
!
!    John Burkardt
!

  

!  integer ( kind = 4 ), parameter :: n = 100
!  integer ( kind = 4 ), parameter :: nc = 51
   ! Need double precision
  integer, parameter :: WP = selected_real_kind(15,307)
  ! Execution status
  integer :: status = 0, ignored_status
  ! Data arrays

  integer ( kind = 4 ), parameter :: n = 256
  integer ( kind = 4 ), parameter :: nc = 129  
  character (len = 256 ) csv_file_name
  real ( kind = 8 ) r8_uniform_01
  integer ( kind = 4 ) i
  real ( kind = 8 ) in(n)
  real ( kind = 8 ) in2(n)
  complex ( kind = 8 ) out(nc)
  
  ! integer ( kind = 4 ) seed
  real*8 samplesec, sec(256),pi,secakhir
  type(DFTI_DESCRIPTOR), POINTER :: hand
  
  ! Data arrays
  real(WP), allocatable :: x_real (:)
  complex(WP), allocatable :: x_cmplx (:)
  
  
  hand => null()

  
  
!  seed = 123456789
  
  csv_file_name = "datainput_mkl.log"
  open(unit=20,file=csv_file_name,status='replace',access ='sequential',form='formatted',action='write')  

  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) 'TEST02'
  write ( 20, '(a)' ) '  Demonstrate mkl on a single vector'
  write ( 20, '(a)' ) '  of real data.'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Transform data to fft coefficients.'
  write ( 20, '(a)' ) '  Backtransform FFT coefficients to recover '
  write ( 20, '(a)' ) '  data.'
  write ( 20, '(a)' ) '  Compare recovered data to original data.'
!
!  Set up the input data, a real vector of length N.
!
!  call r8vec_uniform_01 ( n, seed, in )
!
  samplesec = 0.001d0
  pi = acos(dble(-1.0))
  do i=1,n
      sec(i) = 0.d0 + samplesec*(i-1)
      secakhir = sec(i)
      !amplitudo(i) = sin(2*pi*50*sec(i)) + sin(2*pi*120*sec(i));
      in(i) = sin(2*pi*50*sec(i))
  enddo
 
  
  
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Input Data:'
  write ( 20, '(a)' ) ' '

  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, in(i)
  end do
!
!  Set up a plan, and execute the plan to transform the IN data to
!  the OUT FFT coefficients.
!
  
!  call dfftw_plan_dft_r2c_1d ( plan_forward, n, in, out, FFTW_ESTIMATE )

!  call dfftw_execute ( plan_forward )

  print *,"Example basic_dp_real_dft_1d"
  print *,"Forward-Backward double-precision real-to-complex",               &
          " out-of-place 1D transform"
  print *,"Configuration parameters:"
  print *,"DFTI_PRECISION              = DFTI_DOUBLE"
  print *,"DFTI_FORWARD_DOMAIN         = DFTI_REAL"
  print *,"DFTI_DIMENSION              = 1"
  print '(" DFTI_LENGTHS                = /"I0"/" )', N
  print *,"DFTI_PLACEMENT              =  DFTI_NOT_INPLACE"
  print *,"DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX"

  
  print *,"Create DFTI descriptor"
  status = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_REAL, 1, N)
  if (0 /= status) goto 999

  print *,"Set out-of-place"
  status = DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
  if (0 /= status) goto 999

  print *,"Set CCE storage"
  status = DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE,                   &
                        DFTI_COMPLEX_COMPLEX)
  if (0 /= status) goto 999

  print *,"Commit DFTI descriptor"
  status = DftiCommitDescriptor(hand)
  if (0 /= status) goto 999

  print *,"Allocate data arrays"
  allocate(x_real(N))
  allocate(x_cmplx(N/2+1))

  !print *,"Initialize data for real-to-complex FFT"
  x_real = in
  
  print *,"Compute forward transform"
  status = DftiComputeForward(hand, in, x_cmplx)
  if (0 /= status) goto 999

  
  write (*,*) 'FFT Forward done'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Output FFT Coefficients:'
  write ( 20, '(a)' ) ' '

  do i = 1, nc
    write ( 20, '(2x,i4,2x,2e14.6)' ) i, x_cmplx(i)
  end do
!
!  Set up a plan, and execute the plan to backtransform the
!  complex FFT coefficients in OUT to real data.
!
!  call dfftw_plan_dft_c2r_1d ( plan_backward, n, out, in2, FFTW_ESTIMATE )

!  call dfftw_execute ( plan_backward )
  print *,"Compute backward transform"
  status = DftiComputeBackward(hand, x_cmplx, in2)
  if (0 /= status) goto 999

    !in2 = x_real
  
  write (*,*) 'FFT Backward done'
  write ( 20, '(a)' ) ' '
  write ( 20, '(a)' ) '  Recovered input data divide by N:'
  write ( 20, '(a)' ) ' '

  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, in2(i) / real ( n, kind = 8 )
  end do
!
!  Release the memory associated with the plans.
!
!  call dfftw_destroy_plan ( plan_forward )
!  call dfftw_destroy_plan ( plan_backward )
100 continue
    
  print *,"Release the DFTI descriptor"
  ignored_status = DftiFreeDescriptor(hand)

  
  write ( 20,*) 'difference'
  do i = 1, n
    write ( 20, '(2x,i4,2x,g14.6)' ) i, (in2(i) / real ( n, kind = 8 )) - in(i) 
  end do

  if (allocated(x_real) .or. allocated(x_cmplx)) then
    print *,"Deallocate data arrays"
  endif
  if (allocated(x_real))     deallocate(x_real)
  if (allocated(x_cmplx))    deallocate(x_cmplx)

  stop
  999 print '("  Error, status = ",I0)', status
  goto 100

  
  close(20)
end

Pramudito__Ary · ‎11-10-2017

Problem solved! included PSD, RMS, Full Transform, Half Transform

https://vibrationdata.wordpress.com/2012/10/29/fourier-transform/

http://www.vibrationdata.com/fortran/FOURIER.F

http://www.vibrationdata.com/fortran/FFT.F

Power spectrum fft