How to plot time dependent energy spectrum
 

Hey CFD-Community,

I've done a LES free stream simulation and now i want to plot my data. I have the velocity data of 24 points of my simulation environment for the x, y and z components and i want to plot my data in Matlab. The last plot i need is an energy spectrum plot.

I'm imagining 24 energy spectrum plots for my 24 observer points. But at first i just wanted to a basic energy spectrum plot for one velocity component to get an idea how to do it. I used a method i found here in the forum:

u=fft(VelocityU(:,1));
p=(u.*conj(u))/N^2;

But im not getting any viable data, my energy spectrum has a peak at point 0 and thats it. I have uploaded my result and the origin velocity data. Can you tell me what im doing wrong?

Thank you very much in advance, and have a good day.

First of all, use the log scale but it is also mandatory that you plot the range of frequencies only up to Nyquist cut-off (pi/dt). In a time-frequency analsysis you should adopt a time-windowing technique, why don't you start first by doing a spatial analysis?

Thanks FMDenaro for your reply.

I still have a few questions.
I have implemented the Han-Window and I'm using the log scale now. I'm curious, why should i start with a spatial analysis when I'm not interested in it?(1) Anyway, i did one for my first datarow for my velocity u.
Sidenote: My simulation time is 1 Second with dt=0.0001. Because i want to view a stationary process, i cutout the first 0.2 Seconds, so my data starts with second 0.2001. To fullfil the 2^x law i added 192 zeros to my data (to get to 8192 entries) to use the FFT algorithm. But i think this is a mistake, isnt it? In my time-energy spectrum plot i have small fluctuations for high frequencies, this is due to my addition, isnt it?(2)
Additional, I have a more general question to this transformation. Am I really getting the frequencies, or am i getting the wavenumber here?(3)
My last question for now applies to the spational energy analysis. I have 24 Monitoring points with equal distance to each other. Then my nyquist frequency is just 12?(4)

Some answers were already addressed in similar posts.
However, the window technique should allow you to get some samples of "periodic" temporal signal so that you can do the spectral analysis of each one and then performing a statistical averaging of the computed spectra.


The FFT must be done only for the resolved frequencies, until the Nyquist frequency pi/dt. Use your computational dt to compute.



For a statistically steady flow you have to wait until the numerical transient is finished, only after you collect the time-samples.


The relation between frequency (f) and wavenumber (n) is f=n*2*pi/T. Therefore you need to evaluate the period T.

Thanks again for the respond. I've read a lot of posts in this forum, but sometimes i wasnt able to understand it properly, so i thought i post my specific problem to understand the functinality of the spectrum analysis on the basis of an example i have spend a good amount of time with.

Okay, so my mistake isnt the addition of the zeros in my velocity array, but the neglection of the averaging process. I generally know the functionality of a window function, but im confused of the statistical averaging process. In my first impression i thought you mean a time averaging process by this. So I computed a variable, which is defined as the sum of my fft-transformed velocityvector multiplicated by dt and finally divided with T. T in my case is equal to 0.8192 (N(=8219)*dt(0.0001). But I think this is not what you meant?


I dont get your definition of the nyquist frequency. The nyquist frequency is f_sample/2, and f_sample i can define as 1/Delta_t to sample every datapoint. Am I wrong here? So im getting a nyquist frency of 5.000Hz(Delta_t=0.0001 -> f_sample=10.000Hz). My confusion with the "spatial nyquist frequency" stems from the fact, that if im considering the space and I "hold the time", my data is reduced to the 24 monitor points in my simulation, in contrast to the 10.000 values if I "hold the space" and analyse the time.

Yeah i found thas definition in this forum, too. So my x-axes is wrong, i now have the wavenumber and i need to transform it to the frequencyspectrum.

Thanks a lot for your help so far :) Im making progress.

You don't have to padd the velocity vector with zeros but you have to use the correct dimension. The nyquist frequency (dimensional) is pi/dt and corresponds to the max wavenumber (a number) of your Fourier representation, don't make confusion. In your case you can write pi/dt=n_max*2*pi/T -> n_max=0.5*T/dt is the max wavenumber.



Then, T is not necessarily the total time of your simulation but that characteristic physical period of your flow. If you use the total time of the computation, you are implicitly assuming that the flow repeat yourself every T time. This way you need of N*T times samples to perform the statistical averaging in time.

Okay, I'm really confused right now :confused: I dont get why i would want to know T or the wavenumbers. Im interested in the frequencies. Im doing an fft in the time-domain, so I would say im getting my result in the frequencydomain. I dont need any transformationen to the wavenumberdomain, do I?

And Im still confused about the nyquist frequency defintion. Is my nyquist definition wrong with 5000Hz? [Delta_t=0.0001 -> f_sample=1/Delta_t=10.000Hz -> f_nyquist=f_sample/2).
Because if I use your definiton, im getting a f_nyquist_2=pi/Delta_t=31416Hz.

The Nyquist wavenumber can be defined by the ratio between the max wavelenght T and the smallest one representable on a finite grid of size dt. This is

n_max=T/2dt = N*dt/2*dt=N/2.

The Nyquist frequency is

f_max=n_max*2*pi/T= pi/dt

Thank you very much FMDenaro, your derivation makes perfectly sense. I'm still wondering, why my way of defining the nyqust frequency is wrong? I've tried out both definitions and both lead to the same graph but with a different scaling (see attachments). This is due to my definition of my function, which has an dependency on f_sample:

P1Y =Y_f(1:N/2+1);
PY=(1/(f_sample*N)) * abs(P1Y).^2;
PY(2:end-1)=2*PY(2:end-1);

Does anyone know why my first way of calculating f_nyquist is wrong? I really want to know my error.

As you see, your f_nyquist is actually n_max, that is a number, not a frequency. From your f_nyquist you have to evaluate the cut-off frequency.

Aaah yes now i get it. Sorry, thank you very much :)

Another topic, you said earlier i dont need to add zeros to my velocity vector. I had done this, because i want to get 8192 entries (2^13), so that matlab uses the FFT algorithm and not the DFT one. Im wondering, if i corrupted my data with this approach and now the results are wrong and I maybe shouldn't do the addition off the zeros and instead use the DFT algrotihm?

Im havent seen so many spectrum diagramms until now, so I'm not entirely sure my diagramms are totally correct. Especially the higher frequencies in my diagramms are irritating me, but i guess this is due numerical problems with the simulation.

As reported in mathworks:


"You can interpolate the DFT by zero padding. Zero padding enables you to obtain more accurate amplitude estimates of resolvable signal components. On the other hand, zero padding does not improve the spectral (frequency) resolution of the DFT. The resolution is determined by the number of samples and the sample rate."


https://it.mathworks.com/help/signal...o-padding.html

Again, thank you very much FMDenaro :)

Finally I would like to come back to the last open question in this thread regarding the spatial fft analysis. In my understanding the situation is like this:

I have 24 spatial monitoring Points, so my N_sp=24. So n_max_sp=N_sp/2=12. The corresponding nyquist frequency is f_ny_sp=pi/dx =31.4150Hz because the distance between the points is 0.1m. So to improve my resolution in the space domain i would need to install more monitor points. If im doing everything right, i should be able to see the -5/3 slope in the result diagramm(in the wavenumber domain). Finally, when i want to know my total kinetic energy for every time step E is evaluated like this: E_sp=sqrt(ku^2+ky^2+kz^2). ku,ky and kz are hereby my fft transformed velocity signals. With this i would get for all my 8192 time steps a spatial energy plot, which i can finally sum to the total energy: E_sp_total=Sum_i(E_sp_i*dt)/t_total with i={1,...,8192} which results to the total energy density spectrum, avareged over time. Here, you should be able to see -5/3 slope, too (again, wavenumber domain).

Am i getting something wrong here?

You need to have at least one direction of homogeneity, for example a spanwise periodic direction. Then You collect all the grid points of tour solution for the 1d FFT

Okay, assuming my x direction is homogenous. I still can do a spatial fft with the monitoring points, or am Im wrong? Then my example would be correct if i edit it like this:

I have 24 spatial monitoring Points, so my N_sp=24. So n_max_sp=N_sp/2=12. The corresponding nyquist frequency is f_ny_sp=pi/dx =31.4150Hz because the distance between the points is 0.1m. So to improve my resolution in the space domain i would need to install more monitor points. If im doing everything right, i should be able to see the -5/3 slope in the result diagramm(in the wavenumber domain). Finally, when i want to know my total kinetic energy for every time step E is evaluated like this: E_sp=ku^2. ku is hereby my fft transformed velocity signal. With this i would get for all my 8192 time steps a spatial energy plot, which i can finally sum up to the total energy: E_sp_total=Sum_i(E_sp_i*dt)/t_total with i={1,...,8192} which results to the total energy density spectrum, avareged over time. Here, you should be able to see -5/3 slope, too (again, wavenumber domain).

It makes no sense using the probe points, you have to use the whole nodes of the numerical solution.

The problem is that my grid resolution changes, the further the flow is from the opening the coarser the grid becomes.

If your flow problem has inflow/outflow, you have walls, you have no homogeneity and it makes no sense a spectral analysis. I think you can only do the temporal analysis.

Thank you very much for all your help FMDenaro :)

Perhaps I should have mentioned that I am simulating a free jet into a large volume. But Im gonna stay now with the temporal analysis.

Ok, but did you set periodic conditions in the other two directions? If yes you can do the 1d spectra