Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
26734 Discussions

## Power spectrum fft Beginner
360 Views

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.

6 Replies Valued Contributor II
360 Views

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. New Contributor I
360 Views

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>

``` Beginner
360 Views

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)

!     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 = out*out;  /* 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);
}``` Valued Contributor II
360 Views

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. Beginner
360 Views

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:
!
!
!  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:
!
!
!  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 ``` Beginner
360 Views 