0
$\begingroup$

I am solving a very long problem and one part of it requires me solving an ODE and computing FFT of the resultant data. Essentially I have a differential equation$\frac{dz}{dt}$ for velocity, from which I can derive the acceleration equation. I can then solve the ODE and input the values in the acceleration equation and perform Fourier Transform, or I can take FT of velocity data and multiply it by $\omega$, or I can directly perform FFT on $z$ and multiply the output by $\omega^2$. All these three methods should produce similar spectra.

However, I am getting completely different spectra using these three approaches. I am posting what I think is the relevant part of the code(pertaining to FFT)and the relevant plots. I am not adding complete code to avoid cluttering and possible distraction from the main issue. If someone needs more explanation let me know.

PS: $\beta$ and $\tilde b$ are just two parameters in my ODE. The peak frequency depends inversely on $\tilde b$

Plots: enter image description here enter image description here enter image description here

As requested in the comment section, I am adding the plot of position, velocity, and acceleration in the time domain. enter image description here enter image description here enter image description here Here is the relevant code snippet.

 '''
Fourier Transformation of the acceleration
'''
#https://www.mathworks.com/matlabcentral/answers/15770-scaling-the-fft-and-the-ifft

fft_scaling_factor = 1
spectrum_acc =[]
spectrum_vel =[]
spectrum_pos =[]
velocitylist =[.24]
impact_parameter = [0.01,0.1]

for index,b in enumerate(impact_parameter):
    filter_coeff =[]
    spectrum_a =[]
    spectrum_v =[]
    spectrum_z =[]


    for j in range(len(velocitylist)):
        fft_scaling_factor = 1
        filtered_acc =[]
        filtered_vel =[]
        filtered_pos =[]

        fft_input_acc = acc_normalized[index][j]
        fft_input_vel = velocity_normalized[index][j]
        fft_input_pos = z_list[index][j]

        fft_input_acc = detrend(fft_input_acc,type = 'constant')
        fft_input_vel = detrend(fft_input_vel,type = 'constant')
        fft_input_pos = detrend(fft_input_pos,type = 'constant')



        ft_acc_normalized = rfft(fft_input_acc)*fft_scaling_factor
        ft_vel_normalized = rfft(fft_input_vel)*fft_scaling_factor
        ft_pos_normalized = rfft(fft_input_pos)*fft_scaling_factor

        ft_acc_normalized = np.abs(ft_acc_normalized)**2
        ft_vel_normalized = np.abs([a*b for a,b in zip(ft_vel_normalized,freq_normalized)])**2
        ft_pos_normalized = np.abs([a*b**2 for a,b in zip(ft_pos_normalized,freq_normalized)])**2


        spectrum_a.append(ft_acc_normalized)
        spectrum_v.append(ft_vel_normalized)
        spectrum_z.append(ft_pos_normalized)


    spectrum_acc.append(spectrum_a)
    spectrum_vel.append(spectrum_v)
    spectrum_pos.append(spectrum_z)
$\endgroup$
15
  • $\begingroup$ For me it is not clear where you do your multiplications with $\omega$. Also, the FFT of MATLAB has the 0-frequency at its left, the Nyquist frequency in its middle and 2*Nyquist frequency at the right (which is equivalent to the DC component). This needs to be taken into account when multiplying by $\omega$. $\endgroup$
    – M529
    Dec 9 '19 at 19:09
  • $\begingroup$ I would expect, that if your acceleration spectrum is flat in frequency, then your velocity spectrum would drop as the frequency increases, not rise, and that the position spectrum would drop even faster. Why is this not happening? Or is that your question. It further appears that this "omega arithmetic" is purely the invention of one guy, and is pretty much just the integration rule (i.e., $\mathcal{F}\left \lbrace \frac{d}{dt} x(t) \right \rbrace = \omega \mathcal{F}\left \lbrace x(t) \right \rbrace$), which is hardly new. $\endgroup$
    – TimWescott
    Dec 9 '19 at 19:39
  • $\begingroup$ I deleted my answer to attract a better one and leaving my remaining question here: The description of your math required is correct but I don't quite see where you are multiplying by the frequency and frequency^2 in your code. Also you want to multiply the complex fft output by frequency and then take the PSD of that result; not clear to me that you are doing that? $\endgroup$ Dec 9 '19 at 21:33
  • $\begingroup$ Ok I follow now the multiplication by the frequency axis and your final power computation-- what are the units for frequency in the freq_normalized arrays? $\endgroup$ Dec 9 '19 at 23:16
  • $\begingroup$ @M529 I am using Numpy's rfft so the first frequency term is zero and the last is Nyquist frequency. docs.scipy.org/doc/numpy/reference/generated/… $\endgroup$
    – Prav001
    Dec 10 '19 at 22:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.