Numerical inversion of characteristic functions I have a need to use the FFT in my work and am trying to learn how to use it. I am beginning by attempting to use the FFT to numerically invert the characteristic function of a normal distribution. So I have discretised the integral using the trapezoidal rule (I know, this is a very crude method), converted to a form consistent with the FFT and then run a programme to make the calculations. However when I plot the output, the density function gets thinner and thinner as I increase the number of discretisation steps. I don't know if this is because of my errors or because of the numerical problems associated with the trapezoidal rule. Would someone mind having a look at my working please?
Thanks...
$$
f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-ix\xi}F(\xi)d\xi = \frac{1}{\pi}\int_{0}^{\infty}e^{-ix\xi}F(\xi)d\xi\\
\approx \frac{1}{\pi}\int_0^{\xi_{max}}e^{-ix\xi}F(\xi)d\xi
$$
Let $\xi_j=(j-1)\Delta\xi,\quad j=1,...,N$ 
Take $\xi_N=\xi_{max}$ and $\Delta\xi=\frac{\xi_{max}}{N-1}$ 
Set $F(\xi_j) = F_j$ 
Let $x_k = x_{min} + (k-1)\Delta x,\quad k=1,...,N \quad$ ($\Delta x$ will be determined later)
Set $f(x_k) = f_k$
Discretising integral using trapezoidal rule:
$$
f_k \approx \frac{1}{\pi}\int_0^{\xi_{max}}e^{-ix_j\xi}F(\xi)d\xi \\
= \frac{1}{\pi}\Delta\xi \left(\sum_{j=1}^{N}e^{-ix_k\xi_j}F_j - \frac{1}{2}e^{-ix_k\xi_1}F_1 - \frac{1}{2}e^{-ix_k\xi_N}F_N  \right) \\
= \frac{1}{\pi}\Delta\xi \left(\sum_{j=1}^{N}e^{-i\Delta x \Delta\xi (j-1)(k-1)}e^{-i x_{min}(j-1)\Delta \xi}F_j - \frac{1}{2}F_1 - \frac{1}{2}e^{-ix_k\xi_{max}}F_N  \right) \\
= \frac{1}{\pi}\Delta\xi \left(\sum_{j=1}^{N}e^{-i\frac{2\pi}{N} (j-1)(k-1)}e^{-i x_{min}(j-1)\Delta \xi}F_j - \frac{1}{2}F_1 - \frac{1}{2}e^{-ix_k\xi_{max}}F_N  \right)
$$
where in the last step $\Delta x \Delta{\xi}$ has been set to $\frac{2\pi}{N}$ in order for the sum to be in the form required by FFT. Rearranging gives $\Delta x = \frac{2\pi}{N\Delta \xi} = \frac{2\pi(N-1)}{N\xi_{max}}.$
To centre about the mean $\mu$ set $x_{min} = \mu - \frac{N\Delta x}{2} = \mu - \frac{\pi}{\Delta \xi}.$

Check the normalization (1 or 1/N**0.5 or 1/N) on your FFT solver and make sure you are consistent with that. Also, your end points should be F0 and F(N+1), not F1 and FN.
Here's the solution in Mathematica code:
logCF[k_] := -0.5 k^2
n = 2^10 - 1;
dx = 10/n;

x0 = -dx n/2;
xs = Table[x0 + i dx, {i, 0, n}];
dk = 2 \[Pi]/(n dx);
fs = Table[Exp[logCF[i dk] - I i dk xs[[i]]] If[i == 0 || i == n, 1/2, 1] dk n/\[Pi], {i, 0, n}];
ps = InverseFourier[fs, FourierParameters -> {1, 1}];

Table[{xs[[i]], Re[ps[[i]]]}, {i, n}];
Show[ListPlot[%, Joined -> True, PlotStyle -> {Red}], Plot[PDF[NormalDistribution[0, 1], x], {x, -4, 4}]]

