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) enter image description here enter image description here

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))

        acc_sum =np.add(ft_acc_normalized,acc_sum )
        intensity_list = [power_spectrum_factor * (a ** 2) for a in ft_acc_normalized]
        intensity_sum = np.add(intensity_list,intensity_sum)
    intensity_summed_over_velocity.append(intensity_sum *acceleration_factor**2)

#intensity_summed_over_velocity=  ma.masked_less_equal(intensity_summed_over_velocity,1e-5)

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.ylabel(r'$ a(\tilde t) $', fontsize=14)
    plt.xlabel(r'$ \tilde t $', fontsize=14)
    if screening is False:
        plt.title('a(t) vs time without screening',fontsize=15)
        plt.title('a(t) vs time with screening',fontsize=15)

    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")

    if screening is False:
        plt.title('Single particle spectrum without screening')
        plt.title('Single particle spectrum with screening')


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.

| improve this answer | |
  • $\begingroup$ 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. $\endgroup$ – Prav001 Oct 17 '19 at 20:49
  • $\begingroup$ 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 ? $\endgroup$ – Fat32 Oct 17 '19 at 20:54
  • $\begingroup$ 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. $\endgroup$ – Prav001 Oct 17 '19 at 21:01
  • $\begingroup$ ok so what's is the expected (typical) functional form of this acceleration ? $\endgroup$ – Fat32 Oct 17 '19 at 21:04
  • $\begingroup$ do you have an expression for the dependence of acceleration on $\beta$ abd b 's. $\endgroup$ – Fat32 Oct 17 '19 at 21:08

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.