# The derivative of a gauss function via FFT and IFFT in Python

I have a problem with computing a derivative of a Gauss function using FFT and IFFT from NumPy library. I use the fact that $$$$\frac{d}{dx}f(x) = \frac{1}{\sqrt{2\pi}}\int{ike^{ikx}\mathscr{F}(f(x))(k)dk},$$$$ where $$\mathscr{F}(f(x))(k) = \frac{1}{\sqrt{2\pi}}\int{f(x)e^{-ikx}dx}$$. $$\\$$ My code:

import numpy as np
import scipy as sci
import matplotlib.pyplot as plt

def Gauss(x, mean, sigma):
#GAUSS FUNCTION WITH ITS NORMALIZATION FACTOR
return sci.pi**(-1/4) * sigma**(-1/2) * np.exp(-(x - mean)**2 / (2 * sigma**2))

def Fft(f):
#THAT dx AND np.sqrt(2*sci.pi) ARE SOME NORMALIZATION FACTORS
return np.fft.fftshift(np.abs(np.fft.fft( f ) * dx / np.sqrt(2*sci.pi)))

def Ifft(fftf):
#THAT dx AND np.sqrt(2*sci.pi) ARE SOME NORMALIZATION FACTORS
return np.abs(np.fft.fftshift(np.fft.ifft( fftf )).real) / dx * np.sqrt(2 * sci.pi)

def d_Gauss(x, sigma):
#FINAL DIFFERENCE, LATER DENOTED AS EXACT DERIVATIVE
return (Gauss(x+dx, 0, sigma)-Gauss(x, 0, sigma))/dx

x      = np.linspace(-20, 20, 2**12)
dx     = np.abs(x[1] - x[0])
freq_x = np.fft.fftshift(np.fft.fftfreq(np.shape(x)[0])) / dx * 2 * sci.pi
dk     = np.abs(freq_x[0]-freq_x[1])
sigma  = 3/2
mean   = 0

#COMPUTING EXACT GAUSS DERIVATIVE AND IT'S APPROXIMATION VIA FFT AND IFFT
d_gauss_exact = d_Gauss(x, sigma)
d_gauss       = Gauss(x,mean,sigma)
d_gauss       = Fft(d_gauss)
d_gauss       = 1j*freq_x*d_gauss
d_gauss       = Ifft(d_gauss)

#PLOTTING
plt.figure('Gauss der.')
plt.plot(x, d_gauss,      'b--', label='Gauss der. fft')
plt.plot(x, d_gauss_exact,'r-',  label='Exact Gauss der.')
plt.xlim(-10,10)
plt.xlabel(r'$$x$$')
plt.ylabel(r'$$\frac{d}{dx}g(x)$$')
plt.legend()
plt.show()


What I get is:

The code works fine for FFT of Gauss function, modulated pulse and Lorentz function. That derivative is crucial for me in the next part of the project. I'd be grateful for any hints and help.

In your Ifft function, you are taking the absolute value of the real part. Just take the real part without doing the absolute value.

Here is a correct working Matlab Code for even N:

sigma=3/2;mean=0;N=128;L=20;
dx=L/N;
x=linspace(-L/2,L/2-dx,N);
f=(pi.^(-1/4).*sigma.^(-1/2).*exp(-(x-mean).^2./(2*sigma^2)));
k(1:N/2)=1i/L*2*pi*(0:N/2-1);k(N/2+2:N)=1i/L*2*pi*(-N/2+1:-1);k(N/2+1)=0;
fhat=fft(f);
fhatx=k.*fhat;
fx=ifft(fhatx);
plot(x,f,x,fx);