2

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}.$

Freakalien
  • 229
  • 1
    Are you making some assumption about the symmetry of $F$? How do you justify your symmetry step, right at the beginning? – Nick Peterson Aug 09 '13 at 03:42
  • Hi @NicholasR.Peterson. Sorry for taking so long to get back. f is real (its a density function) so F(-y) = F*(y) which means that F is even in its real part and odd in its imaginary part. – Freakalien Aug 10 '13 at 00:55
  • Hi @NicholasR.Peterson. I have tried computing f with a two-sided integral but the phenomenon remains (although to only about half the extent). – Freakalien Aug 10 '13 at 03:45

2 Answers2

1

Please have a look the section on the numerical inversion of characteristic functions in the dissertation found here at 6.3.1.3:

http://hdl.handle.net/10539/9273

The algorithm is outlined and there is Matlab code

It maybe helpful.

Renay Singh
  • 11
1

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}]]
Steve Schulist
  • 107
  • 7