Signal response amplitude depends on the time interval in simulation

I have already tried to look for an answer, but I do not find existing answers satisfactory.

I am interested in the absolute value of the response function of a damped oscillator (or any time series). However, when I change my simulation interval, the response amplitude which I obtain with a Fourier transform, also changes. In my case I am interested in the absolute value because I need to relate it to physical parameters in my system, and I need to obtain something which is fixed.

See attached image of time series and the respective Fourier transforms for different time intervals.

Here is also the code if you want to simulate:

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp
import scipy.fft as fft

tsteps = 10000
tmax=1000
t = np.linspace(0, tmax, tsteps)
init_condition=[1, 0]

#solving harm oscillator
def harm_osc():
def dndt(t, V):
x=V[0]
xdot=V[1]
gamma, omega0 = 1, 2
xdot2=-gamma*xdot-omega0**2*x
return [xdot, xdot2]

sol = solve_ivp(dndt, [0, tmax], init_condition, method='LSODA', t_eval=t)
return sol

sol=harm_osc() #solve the equation
plt.plot(sol.t, sol.y[0])
plt.xlabel("time")
plt.ylabel("amplitude")
plt.show()

normalization=2/tsteps
sol_fft=np.abs(fft.fft(sol.y[0])*normalization) #normalize the amplitude with timesteps
sol_fft=sol_fft[0:tsteps//2] #retaining only positive frequency values
freq = fft.fftfreq(tsteps, 1/tsteps)
freq=freq[0:tsteps//2]#retaining only positive frequency values

plt.plot(freq, sol_fft)
plt.xlim(1e-1,1e1)
plt.xscale("log")
plt.xlabel("frequency")
plt.ylabel("amplitude")
plt.show()



Is there a way to fix it? I am interested directly in the amplitudes of my oscillator, rather than power. In this case, that would be a displacement of the oscillator.

• Why did you choose this normalization? Seems wrong for an impulse response. Commented Apr 18 at 16:40
• @Hilmar I read about the normalization here: towardsdatascience.com/…. Obviously it is not what I want, but right now I do not know any better. That is basically what my question asks - how to find the correct normalization Commented Apr 18 at 17:07
• If you use the DFT to convert an impulse to a transfer function and vice versa, don't normalize at all. Just use the "text book" DFT convention without any scaling. Commented Apr 18 at 18:05
• @Hilmar Sorry, but that is exactly the problem - if I do not normalize, the amplitude depends on my time interval. The longer time I take, the smaller the amplitude. Commented Apr 18 at 18:54

Reading your question and comments, and testing your code, my understanding of your problem is that you do not understand why you need to normalise the transformed signal so that the amplitudes in the frequency domain equal those in the time domain, and that the Python code you shared solved the problem, but you don't understand why.

With that in mind, let me try to demonstrate why that particular normalisation you used in your Python code is required by working through the mathematics of the Fourier transform using a simple example.

Consider a signal that is $$x[n]=A\cos(\omega n)$$. Referring to the scipy.fft documentation, you can see the definition of the 1D discrete Fourier transform is that it multiplies your signal by a complex exponential: $$y[k]=\sum^{N-1}_{n=0} x[n]e^{-j\omega_kn}=A\sum^{N-1}_{n=0}\cos\left(\omega n\right)\left(\cos\left(\omega_kn\right)-j\sin\left(\omega_kn\right)\right)$$ Where I have used Euler's formula and: $$\omega_k=\frac{2\pi k}{N}$$ Then, using product-to-sum trigonometric identities gives: $$y[k]=\frac{A}{2}\sum^{N-1}_{n=0}\left[\cos\left((\omega-\omega_k)n\right)+\cos\left((\omega+\omega_k)n\right)\right]-j\left[\sin\left((\omega-\omega_k)n\right)-\sin\left((\omega+\omega_k)n\right)\right]$$ And so when $$\omega=\omega_k$$, the following results: $$y[k]=\frac{A}{2}\sum^{N-1}_{n=0}\left[1+\cos\left(2\omega_k n\right)\right]-j\left[0-\sin\left(2\omega_k n\right)\right]\\=A\frac{N}{2}+\frac{A}{2}\sum^{N-1}_{n=0}\cos\left(2\omega_k n\right)-j\sin\left(2\omega_k n\right)$$ Now, it is important to recognise that the frequencies ($$\omega_k$$) used in the Fourier transform cause the summation on the right to equal zero, and therefore: $$y[k]=A\frac{N}{2}$$ Hence, to get the amplitude, $$A$$, of non-zero frequencies in the Fourier transform, you must divide by $$\frac{N}{2}$$, which is exactly what your Python code normalisation is. However, note that if you work through the same process for $$\omega_k=0$$, you'll get: $$y[0]=\sum^{N-1}_{n=0} x[n]e^{-j\omega_kn}=\sum^{N-1}_{n=0} Ae^{0}=AN$$ Therefore, the normalisation for $$\omega_k=0$$ instead involves dividing by $$N$$.