I am trying to use the result that $p(omega)^n$ = Gamma(n, beta) where $p(omega)$ is the fourier transform of an $Exp(beta)$ distribution. My code does well for small n but as n increases the gamma starts to move across the domain faster. Here I have just set beta=1 for simplicity. I can not find the error in my code.

Note the mask is just to remove the horizontal Line artefact that often occurs after numerical fourier transforms.

I have pasted the code below.

`import numpy as np
import matplotlib.pyplot as plt
#DETERMINISTIC CASE 

n = 15   # Scaling factor for the Fourier transform
N = 1024     # Number of points
T = 0.01     # Sampling interval
shift_point = 100 # Desired center of the plot

# Discretize the PDF
x_values = np.linspace(xm, 100, N, endpoint=False)
pdf_values = expon.pdf(x_values, scale=1)

# Numerical Fourier Transform of the Pareto PDF
P_omega = np.fft.fft(pdf_values)
# Scale and account for deterministic number of sources 
exp_nP = P_omega**(n)

# Compute the full frequency array for IFT
full_exp_nP = np.concatenate([exp_nP, np.conj(exp_nP[::-1])])

# Inverse Fourier Transform
ift_result = np.fft.ifft(full_exp_nP)
ift_result=ift_result[::2]
# Plotting

#Remove all small values
mask = ift_result > 0.01
ift_result = ift_result[mask]
x = np.fft.fftfreq(N, T)[mask]


plt.figure(figsize=(10, 5))
plt.plot(x, -np.abs(ift_result)/np.trapz(x, np.abs(ift_result)), label = "fourier then inverse fourier ")
plt.title('Inverse Fourier Transform of exp(n * P(w)) Centered at 100 for Pareto Distribution')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.grid(True)
##TEST THE EXPONENTIAL DISTRIBUTION AS AN ALTENRATIVE TO THE PARETO. CHECKING HOW THE SUM OF EXPONENTIALS IS DISTRIBUTED. 

#GAMMA RV
x_gamma = np.arange(0, 50, 0.01)
y_gamma = gamma.pdf(x_gamma, n, scale=1)
plt.plot(x_gamma, y_gamma, label = "Gamma(n, beta)")
plt.legend()
plt.show()`