FFT spectra for fundamental frequency from timeseries data (1D REAL to COMPLEX ) -> (REAL Spectra)

[Pramudito__Ary]

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.

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.

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.

1. 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?

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 ?

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 ??

 !    example data...

      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 get power spectra??? conjugate ??



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


end

!The fftw C examples spectra...how to doit in fotran?

! 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);
}

[MGRAV]

I am not completely sure to answer your questions.

First, you need to know that it's time to change form fftw2 to fftw3, like say one the website it's a 20-year-old library. As you can imagine, performance for new hardware is poor.

fftw3 example in Fortran http://www.fftw.org/fftw3_doc/Fortran-Examples.html

If you need to do multiple times the same fft, with only new value each time (but the same size). You can create a plane with another option and reuse it:

http://www.fftw.org/fftw3_doc/Planner-Flags.html

When you use MKL with the fftw3 interface you need to include fftw3_mkl.h, at least in C.

I seems to me that the planner flag is not important anymore, mkl don't use it. However, I like to have the good one that I can easily switch from mkl to fftw3 without to rewrite the full code.

About computing power, for sure you can multiply by the conjugate, or just take squared norm, R^2+I^2, with R and I, respectively the real and complex part.

PS: Matlab gives I always full complex fft, that mean if the input is real, the result is symmetric. However, in fftw and mkl are smart and compute only the needed part. You finish with only with half of the array.  

PS: some size of arrays are more efficient than others, I let you check in the documentation.