1
$\begingroup$

I have a problem with computing a derivative of a Gauss function using FFT and IFFT from NumPy library. I use the fact that $$ \begin{equation} \frac{d}{dx}f(x) = \frac{1}{\sqrt{2\pi}}\int{ike^{ikx}\mathscr{F}(f(x))(k)dk}, \end{equation} $$ 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:

enter image description here

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.

$\endgroup$
1
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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

enter image description here

$\endgroup$

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.