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

Hi Lourival, If i understood you, you want to calcualte the spectrum of time-series 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

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

It's me again, the FFT(U) returns a result like z=a-ib 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

Dear Lourival, YES, right, if z=x+iy, the conjugate of z is x-iy. 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

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 = 2E-5 [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...

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

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

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 m-file. 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

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

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.

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

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

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

Dear Lourival,

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

Zack

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) sub-segments,

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

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 sub-segments (50% overlaped) For example total data of 4096, i.e. N = 3, the first sub-segment 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/(N-1)) ? On which k is the index of the full data (from 1 - 4096)??

Thanks a lot for your help and attention...

Best regards

Lourival

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 sub-segment: 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

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 N-1, we take segments, possibly overlapping, of length L with starting points of D unity apart.

X1(j) = X(j), j=0..L-1 X2(j) = X(j+D), j=0..L-1 X3(j) = X(j+(K-1)D), j=0..L-1 .... 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..L-1 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 L-1.

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 L-1 Am I not modifing the average spectrum??.

Sorry bother you. And don't worry about english, mine isn't good enought

Thanks

Lourival

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 non-square 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