
[Sponsors] 
May 16, 2007, 12:58 
Fast Fourier Transform

#1 
Guest
Posts: n/a

Hi everyone, I would like to analyse the turbulent spectrum of the flow over a bluff body. As I'm simulating using FLUENT all the data is on the Spatial space and I need to use FFT to transform the data to the Fourier space.
In order to compare my simulation with other results I got the velocity (U, V, W) on the same point that the paper presented. Also I have the code to compute the FFT. So my question is about how to do it. To obtain the spectrum, should I use the FFT of the V velocity and from the results get the square ? For example: FFT(V) = R+iIm The one dimensional streamwise velocity spectrum is R^2+Im^2 plotted in a log graph? What about the frequency on the X axis? How do I calculate? Does any one know about it? Thanks Lourival 

May 17, 2007, 09:35 
Re: Fast Fourier Transform

#2 
Guest
Posts: n/a

Hi Lourival, If i understood you, you want to calcualte the spectrum of timeseries of the velocity. Spectrum of U for example can be computed like this spectrum of U=FFT(U)*conjugate of FFT(U), for V just replace U. Regards, Zack


May 17, 2007, 10:37 
Re: Fast Fourier Transform

#3 
Guest
Posts: n/a

Hi Zack, you got the problem, that is right I want to calculate the spectrum of time series of the velocity.
Sorry about my small knowledge about Fourier, I'm trying to get some more insigths about it, my doubts is about the conjugate. The conjugate of FFT(U) is the imaginary part? So the spectrum of U is the (Real part of FFT(U))*(Imaginary part of FFT(U)) ? Another doubte, when I calculate the FFT(U) I have a simetry on the half way, for example 512 on the 1024 data. Should I consider this in my calculus? Regards Lourival 

May 17, 2007, 12:04 
Re: Fast Fourier Transform

#4 
Guest
Posts: n/a

It's me again, the FFT(U) returns a result like z=aib and the conjugate is z'= a+ib, so the spectrum is = a^2 + b^2, OK??
Also, as I can see, once we calculate the FFT(U) we have a plane of simetry, as I said on the last msg. Consequently the spectrum is related only to the first part (or second). To determine the frequency on the abscissa axis it is necessay to calculate w = arctan(b/a) for each point of the spectrum? Another question is concerned about a filter, is necessary to filter the results? Is necessary to filter the spectrum velocity? Thanks for your attention Best regards Lourival 

May 17, 2007, 13:41 
Re: Fast Fourier Transform

#5 
Guest
Posts: n/a

Dear Lourival, YES, right, if z=x+iy, the conjugate of z is xiy. About the frequaency, i suppose you want to calculate the Strouhal number since you are simulating flow around a bluff body. Am I right?. OK. the frequency depends on sampling rate in a dimensional sence(the rate that write the data to the file per second), frequency is 1/T. be aware to calculate different frequencies for only half of the sepctrum as you said. for filtration, as i know it depends. Regards, Zack


May 17, 2007, 14:25 
Re: Fast Fourier Transform

#6 
Guest
Posts: n/a

Dear Zack... again you are right. As I'm simulating the flow around bluff body I would like to have the Strouhal number.
The values of the velocity were sampled at each time step, so the sampling rate were the same of the time step, in this case ts = 2E5 [s]. Unfortunately I didn't get the idea concerning the frequency. Tell me if I am wrong. The frequency on the abscissa axis (log scale) is 1*1/ts, 2*1/ts, 3*1/ts... and so on??? Thanks again... 

May 17, 2007, 14:55 
Re: Fast Fourier Transform

#7 
Guest
Posts: n/a

Dear Lourival, Yes, you are right. I will give example from Matlab for 512 samples. Yfft = fft(y,512); this fft of the samples Power_spectrum = Yfft.* conj(Yfft) / 512; Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis: f = sample rate*(0:256)/512; from 0 to the half of the samples. Strouhal number =f*length scale/velocity scale. Hope everything is clear now. Regards Zack


May 17, 2007, 16:07 
Re: Fast Fourier Transform

#8 
Guest
Posts: n/a

Dear Zack, I believe that I got the idea. I just didn't understand why Power_spectrum = Yfft.* conj(Yfft)is divided by 512.
So the frequency is f = sample rate*(from 0 to 256)/512 and for each frequency, f , should I calculate the Strouhal number? And then compare each one with the Theoretical Strouhal number? Thanks a lot Lourival 

May 17, 2007, 16:39 
Re: Fast Fourier Transform

#9 
Guest
Posts: n/a

Hi Lourival, First of all, you can calculate for any amount, say to all data N, when you calculate f should be from (0 to N/2)/N. when you done that, plot power spectrum (from 1 to (N/2+1))Vs f (or St), f or St you want is the value when power spectrum is the highest. see here a part from my mfile. UFFT=fft(U,N); spectrum=UFFT.* conj(UFFT)/N; f = sample rate*(0:N/2)/N; St=f*(cylinder diameter)/(U mean at inlet); figure(1) plot(St,spectrum(1:N/2+1)) title('STROUHAL NUMBER') xlabel('Strouhal Number') ylabel('Power Spectrum')
Regards, Zack 

May 18, 2007, 16:49 
Re: Fast Fourier Transform

#10 
Guest
Posts: n/a

Hi Zack, thanks a lot for your help, I did exactly as you told me, I believe that I have a small data on the velocity and this results in a lot of noise and it is not possible to identify the real Strouhal number... I have 4000 data of velocity (I filled all the 4096 data, enought to be proportional to 2, with 0). What do you think about it??
Thanks for your help .. Lourival 

May 20, 2007, 15:58 
Re: Fast Fourier Transform

#11 
Guest
Posts: n/a

Dear Lourival, What is your Reynolds number?, if is high enough to be turbulent, then it is normal to see noise in the results. also you can plot the velocity and see if it is periodic laminar or there are fluctuations. in case of laminar vortex shedding you should get one peak.


May 21, 2007, 11:10 
Re: Fast Fourier Transform

#12 
Guest
Posts: n/a

Hi Zack, well my Reynolds number is on the transition zone, ReD = 3900. I believe that there are noises, specially because on the velocity profile there is fluctuations around a average profile and it is periodic but probably not enough. As I said before I took only 4000 samples of the velocity, and I'm starting to believe that at least the double is the minimum... What do you think?
Best Regards Lourival 

May 21, 2007, 18:27 
Re: Fast Fourier Transform

#13 
Guest
Posts: n/a

Dear Lourival
I see your flow is turbulent. Actually for 2D you can see 2 peaks at Re=>250, so with inceasing Re will get more peaks in strouhal number indicating to the fluctuation, so the shedding has no unique frequancy. Regards, Zack 

May 22, 2007, 07:43 
Re: Fast Fourier Transform

#14 
Guest
Posts: n/a

Thanks Zack, Actually I agree with you, but I belive that there is a peak indicating the Strouhal number, as many others, but the Strouhal should be the expressive one. And I don't get this...
Best Regards Lourival 

May 22, 2007, 09:49 
Re: Fast Fourier Transform

#15 
Guest
Posts: n/a

Dear Lourival,
OK. could please send me your data I will try to see it. thesis_ak at yahoo dot com Zack 

May 29, 2007, 07:12 
Re: Fast Fourier Transform

#16 
Guest
Posts: n/a

Hi,
To reduce the noise level in your spectrum (at the expense of spectral resolution), use Welch's method: 1) divide your data in N (possibly 50 % overlapped) subsegments, 2) for each segment compute the power spectrum using FFT over this segment: S_i(f)=FFT_{u_i(t)}(f)*conjugate[FFT_{u_i(t)}(f)], 3) average the N spectra : S(f)=(1/N)*Sigma_{i=1}^N S_i(f) To improve the results, you should multiply your data by a window (whatever your like : Hamming, Hanning, Welch, Parzen...) prior to compute the FFT in stage 2 If spectrum levels are of importance, you should ensure that your spectrum is normalized as in the experiments. Most probably the normalization is such as the integral from 0 to half the sampling rate of S(f) df is equal to the variance of your data (Parseval theorem). Also note that for this kind of spectrum with strong peaks, it is often interesting to use more advanced spectrum estimators than FFT ones, especially if you are a little short on data. You could for instance try Burg's method (maximum entropy). Details and C or Fortran programs can be found in numerical recipes online. Hope this helps. Lionel 

May 30, 2007, 14:12 
Re: Fast Fourier Transform

#17 
Guest
Posts: n/a

Dear Lionel, I believe that I got the idea of reducing noise on the FFT. Please check if it is right.
1) Divide my data on N subsegments (50% overlaped) For example total data of 4096, i.e. N = 3, the first subsegment goes from 1  2048, the second goes from 1024  3072, the last from 2048  4096. 2) After this sub division calculate the FFT of each sub segment and multiply by its conjugate. 3) Average the N spectra. Well if it is right I got the idea, my doubts are concerning the window, should I multiplicate my data on the sub segment by a function that is validy only on that range? For example, for the first sub segment the function is zero from 2049  4096? Another doubt concerning the FFT of the subsegment, should I calculate the FFT for the first sub segment, for instance, from 1 to 4096 or only for 1  2048? On the averaging I didn't get the meaning of Sigma_{i=1} on your answer. For the window function, the Hann function is w(k) = 0.5*(1 cos(2*k*pi/(N1)) ? On which k is the index of the full data (from 1  4096)?? Thanks a lot for your help and attention... Best regards Lourival 

May 31, 2007, 06:14 
Re: Fast Fourier Transform

#18 
Guest
Posts: n/a

Hi,
There are two possibilities: 1) Use a FFT of length equal to the length of the sub segment: using your example, 2048 => 1025 independent Fourier modes with frequencies ranging from 0 to f_sampling/2. The window should have the same length than the segment, therefore N=2048 in your formula. 2) Use a FFT of length equal to twice (or any other integer multiple) the length of the subsegment: using your example, the data d(i) that have to be Fourier transformed are of length 4096 and are built such as: i=1 > 2048: d(i)=w(i)*s_j(i) with w the window function (N=2048) and s_j the samples over the jth segment (i.e. s_1(1)=s(1), s_1(2)=s(2), ... , s_1(2048)=s(2048), s_2(1)=s(1025), s_2(2)=s(1026), ... , s_2(2048)=s(3072), s_3(1)=s(2049), s_3(2)=s(2050), ... , s_3(2048)=s(4096)) i=2049 > 4096: d(i)=0.0 (known as zero padding). You will obtain 2049 Fourier modes with frequencies ranging from 0 to f_sampling/2. However only half of them, the even ones, are independent (and identical to the ones obtained when using method 1), the other being related to a kind of spectral interpolation. Note that because of the windowing, the power contained in the spectra is altered due to the power of the window. It should be corrected either:  by multiplying each jth spectrum by a constant computed so as to ensure that the Parseval equality is met (integral over spectrum related to the jth segment=variance over this segment) or  by dividing the resulting jth spectrum by W=N*sum_(i=1>N) w(i)^2. The first approach is safer since it is free of difficulty that could arise from the normalization used in the FFT program. Lastly, in my previous post, sigma_{i=1} was simply meaning sum from i=1. Hope this helps (sorry for the bad english) Lionel 

May 31, 2007, 18:03 
Re: Fast Fourier Transform

#19 
Guest
Posts: n/a

Dear Lionel, I got the paper of Welch that it is described the method. Here I resume the method.
"Supose X(j), j=0 to N1, we take segments, possibly overlapping, of length L with starting points of D unity apart. X1(j) = X(j), j=0..L1 X2(j) = X(j+D), j=0..L1 X3(j) = X(j+(K1)D), j=0..L1 .... Xk(j) that cover the entire record For each segment of length L we calculate the modified periodogram. That is select a data window W(j), j=0..L1 and form the sequences X1(j)W(j), ..., Xk(j)W(j). We take de FFT of each segment: A1(n), ..., Ak(n). Finally we obtain the K modified periodogram Ik(fn)=L/U*Abs(Ak(n))^2, k=1,...,K. where: fn = n/L, n=0,...,L/2 And U = 1/L*Sum(W(j)^2) that j=0 to L1. So that the spectral estimate is the average of these periodogram P(fn) = 1/K*Sum(Ik(fn)) that k=1 to K." As I understood correctly the average is over the half part of the record for each segment and is it different from that presented by you? If I take for simplicity W(j) = 1 for j=0 to L1 Am I not modifing the average spectrum??. Sorry bother you. And don't worry about english, mine isn't good enought Thanks Lourival 

June 1, 2007, 08:55 
Re: Fast Fourier Transform

#20 
Guest
Posts: n/a

Hi,
The Welch method, as summarized in your post, seems to me similar to the one I was describing: K segments (possibly overlapping) of length L resulting in an average over K spectra of length L/2+1 (frequencies from 0 to f_sampling/2 or reduced frequency f/(f_sampling/2) from 0 to 1, as in Welch paper). For the sake of simplicity you can use a square window (W(j)=1). However it implies that for each frequency the spectrum you obtain is equal to the convolution of the underlying spectrum with the Fourier transform of the square function (0 outside of the segment and 1 inside) that is equal to sinc(f)=sin(f)/f, a slowly decaying function of f. Consequently, part of the energy of one single Fourier mode will be spread over the adjacent frequencies. This could be a serious problem for spectra with strong peaks as the one you want to analyze. The purpose of a window is to use a function with a Fourier transform that has a faster decay than the sinc function so as to reduce as much as possible this smearing. Moreover, since part of the data has been damped near the beginning and the end of the segment when using a nonsquare window, overlapping can be used to partly recover the influence of the damped samples since they will be used in the computation of more than one spectra. Regards, Lionel 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Fast Fourier Transform  alastormoody11  STARCCM+  0  January 20, 2011 23:49 
Poisson equation fourier transform before discretization  yohey44  Main CFD Forum  0  November 9, 2010 13:00 
Fast Fourier Transform of CL signal  Sachin Paramane  Main CFD Forum  0  June 11, 2007 07:51 
Fast Fourier Transform  Lourival  FLUENT  0  May 16, 2007 14:02 
Fast Fourier Transform  Sara  Main CFD Forum  2  February 14, 2006 08:51 