# Getting non smooth spectrum after performing fft on numerically obtained acceleration data

I am trying to calculate the spectrum of Bremmstrahlung, which involves calculating the Fourier transformed acceleration. I am solving a non-linear ODE to numerically calculate the acceleration in the time domain. After taking the Fourier transformation using Numpy's fft, the resultant spectrum looks highly non-smooth and "non-physical" . I cannot paste the entire code so I am posting what I think is relevant snippet. Can someone point out what I am doing wrong?

Note: My acceleration is a function of two variables (beta, and b, the impact parameter), and I want to plot it for the different b, and in order to factor out the beta term I am just summing over all the values of acceleration for the different beta, for a given impact parameter b.

Also my spectrum is Fourier transformed acceleration square times a constant factor (Larmor's formula)

'''
Fourier Transformation of the acceleration
'''
N = 2**8
sampling_frequency = 10
a_w_normalized = [[] for a in range(len(impact_parameter))]
a_w =[[] for a in range(len(impact_parameter))]
intensity_normalized = [[] for a in range(len(impact_parameter))]
acc_summed_over_velocity = []
intensity_summed_over_velocity =[]

for index,b in enumerate(impact_parameter):
acc_sum =[0 for x in range(N)]
intensity_sum = [0 for x in range(N)]
for j in range(len(velocity_z_component)):
#window_kaiser = signal.kaiser(N, 15)
#window_hann = signal.hann(N,sym=True)
window =1
fft_input = acc_normalized[index][j]*window
ft_acc_normalized = np.abs(np.fft.fft(fft_input,norm=None))

intensity_list = [power_spectrum_factor * (a ** 2) for a in ft_acc_normalized]
intensity_normalized[index].append(intensity_list)
a_w_normalized[index].append(ft_acc_normalized)
acc_summed_over_velocity.append(acc_sum)
intensity_summed_over_velocity.append(intensity_sum *acceleration_factor**2)

'''
Plotting acceleration for selected values of impact paramters in time and frequency domain
'''

if plot is True:
plt.figure(figsize=(12, 8))
for i in range(len(impact_parameter)):
plt.plot(t,acc_timedomain_summed_over_velocity[i], label='b={:.3f}'.format(impact_parameter[i]), )
plt.legend()
plt.ylabel(r'$$a(\tilde t)$$', fontsize=14)
plt.xlabel(r'$$\tilde t$$', fontsize=14)
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_smart_bounds(True)
ax.spines['bottom'].set_smart_bounds(True)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
if screening is False:
plt.title('a(t) vs time without screening',fontsize=15)
plt.savefig('./Plots/acc_without_screening.png')
else:
plt.title('a(t) vs time with screening',fontsize=15)
plt.savefig('./Plots/acc_with_screening.png')

'''
Fourier transformed intensity for different impact paramater
'''

plt.figure(figsize=(12, 8))
freq_normalized = np.fft.fftfreq(N)*(2*np.pi*sampling_frequency)
for index,b in enumerate(plasma.impact_parameter):
plt.plot(np.abs(freq_normalized), intensity_summed_over_velocity[index], label=r'$$\tilde b={:.2f}$$'.format(b), )
plt.xlabel(r'$$\tilde \omega$$', fontsize=14)
plt.ylabel(r'$$I_{\omega}$$', fontsize=16)
plt.legend(loc="lower left")
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

if screening is False:
plt.title('Single particle spectrum without screening')
plt.savefig('./Plots/spectrum_no_screening.png')
else:
plt.title('Single particle spectrum with screening')
plt.savefig('./Plots/spectrum_debye_screening.png')

plt.show()


Most of the time, a non-parametric PSD estimation based on FFT alone (called a Periodogram) will provide a random looking spectrum of the numerical data.

To reduce those random variations, or to get a smoother looking estimate, you can do one of the two things:

• Use Periodogram averaging such as Welch's method.
• Use a model based parametric PSD estimate instead.

Note that the first approach will trade-off spectral resolution for reduced variance, whereas the second method will be high resolution but requires a consistent knowledge of the associated physical model of the process.

• Sorry if my comments doen't make much sense, this is the first time I am reading about PSD, and Welch's method. Are you saying that I should calculate my spectrum using PSD instead of squaring Fourier domain acceleration and then multiplying by relevant coefficient? I am interested in calculating Larmor's power which in frequency domain is given by 8*pi*e^2 a(w)^2/3c^3. – Prav001 Oct 17 at 20:49
• I don't know about Larmor's Power $8 \pi e^2 a(w)^{2/3} c^3$... Is $a(w)$ the Fourier transform of acceleration ? – Fat32 Oct 17 at 20:54
• Yes, $a(\omega)$ is the Fourier transformation of acceleration in the time domain. So essentially I am performing FFT on my acceleration data which I show in the question. Now, my acceleration depends on two parameters, b and $\beta$. What I show in my acceleration plot is acceeleration summed over all values of $\beta$ for different b's. – Prav001 Oct 17 at 21:01
• ok so what's is the expected (typical) functional form of this acceleration ? – Fat32 Oct 17 at 21:04
• do you have an expression for the dependence of acceleration on $\beta$ abd b 's. – Fat32 Oct 17 at 21:08