- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page