\chapter{Numerical Integral Transforms}
\label{app:transforms}

In this appendix I want to compile some useful information about the numerical
computation of integral transforms, especially Fourier and Hilbert transforms.
The way we will tackle the computation of these transforms here is by
representing the function to be transformed as a linear combination of so
called B-splines.  The computation of the integral transform is then reduced
to the computation of the integral transform of the B-splines which can be
done analytically.  In this manner we obtain an accurate representation of the
integral transform and in the case of Fourier transforms we will not have to
worry about aliasing.   This approach will also allow us to treat properly any
discontinuities that might be present in our functions.  

Before we start on the computation of transforms let us briefly discuss the
representation of function in terms of B-splines.

\section{B-splines}
B-splines are certain types of polynomial functions that can be used as basis
functions for polynomial splines of functions.  The way we usually think about
splines is to represent our function as a polynomial over all the intervals
that we consider and then the parameters of those polynomial are found by
matching the function value and the derivatives at the interval boundaries,
i.e. at the data points.  With the B-splines, which we here use only for
regular grids, the parameters of the polynomial function being used are fixed
and our function is represented as a linear combination of translations of the
basis function.  I.e.~if we have a function $g(x)$ defined in the interval
$[x_1,x_N]$ with $x_{i} = x_{1}+(i-1)\Delta$, then we represent $g(x)$ within
the interval as:
\begin{equation}
  g(x) \simeq \sum_{i=1}^{N}\tilde{g}_{i}
  \phi\left(\frac{x-x_{i}}{\Delta}\right)
\end{equation}
with the requirement that the following condition holds:
\begin{equation}
  g(x_j) = g_j = \sum_{i=1}^{N}\tilde{g}_{i}
  \phi\left(\frac{x_{j}-x_{i}}{\Delta}\right)
\end{equation}
The only types of B-splines that we will consider here are linear and cubic
B-splines. The normalized linear B-spline is given by the following equation
\begin{alignat}{2}
\phi_1(x) &= 0    &\qquad&\mbox{if } x \leq -1
\nonumber \\
&= 1+x  &\qquad&\mbox{if } -1 \leq x \leq 0  
\nonumber \\
&= 1-x  &\qquad&\mbox{if } 0 \leq x \leq 1  
\nonumber \\
&= 0   &\qquad&\mbox{if } 1 \leq x  
\end{alignat}
The points $\{-1,0,1\}$ are called the knot points of the spline.  The
normalized cubic B-spline on the knots $\{-2,-1,0,1,2\}$ is given by
the following equation
\begin{alignat}{2}
\phi_3(x) &= 0   &\qquad& \mbox{if } x \leq -2 
\nonumber \\
&= {1 \over 6}(2+x)^3 &\qquad& \mbox{if } -2 \leq  x \leq -1  
\nonumber \\
&= {2 \over 3}-x^2-{1 \over 2}x^3 &\qquad& \mbox{if } -1 \leq x \leq 0  
\nonumber \\
&= {2 \over 3}-x^2+{1 \over 2}x^3 &\qquad& \mbox{if } 0 \leq x \leq 1  
\nonumber \\
&= {1 \over 6}(2-x)^3 &\qquad& \mbox{if } 1 \leq x \leq 2  
\nonumber \\
&= 0   &\qquad& \mbox{if } 2 \leq x   
\end{alignat}
We should note here that there is a simple way of generating these B-splines
from the basis function:
\begin{alignat}{2}
\phi_0(x) &= 0    &\qquad&\mbox{if } x \leq -\frac{1}{2}
\nonumber \\
&= 1  &\qquad&\mbox{if } -\frac{1}{2} \leq x \leq \frac{1}{2}  
\nonumber \\
&= 0  &\qquad&\mbox{if } \frac{1}{2} \leq x.  
\end{alignat}
Namely we have:
\begin{equation}
\phi_{1}(x) = \int_{-\infty}^{\infty}\phi_{0}(y)\phi_{0}(x-y)dy
\qquad\mbox{and}\qquad
\phi_{3}(x) = \int_{-\infty}^{\infty}\phi_{1}(y)\phi_{1}(x-y)dy
\end{equation}
and in principle we could generate higher order B-splines in this manner.
\begin{figure}[t]
\centering
\begin{tabular}{c@{\;}c}
\mbox{\epsfig{file=Figures-App/BsLin.ps,angle=-90,width=\halfpage}} &
\mbox{\epsfig{file=Figures-App/BsCub.ps,angle=-90,width=\halfpage}} \\
\end{tabular}
\caption{Examples of B-splines and the overlapping of the splines between
  adjacent mesh points. On left is the linear B-spline and on right is the
  cubic one.}
\label{fig:B-splines}
\end{figure}   

We see clearly from figure \ref{fig:B-splines} that for the linear
B-spline the solutions for the $\tilde{g}_{i}$'s is simply
$\tilde{g}_{i} = g_{i}$ but for the cubic B-spline we need to solve a
set of linear equations.  Namely we have:
\begin{alignat}{2}
\frac{1}{6}\left(4\tg_{1}+\tg_{2}\right) 
&= g_1  & & 
\nonumber \\
\frac{1}{6}\left(\tg_{i-1}+4\tg_1+\tg_{i+1}\right) 
&= g_i  &\quad&i = 2,\ldots,N-1 
\nonumber \\
\frac{1}{6}\left(\tg_{N-1}+4\tg_{N}\right) 
&= g_N.  & & 
\end{alignat}
In this manner we can of course only match the function values but with the
cubic B-spline we should also be able to match the derivatives of the function
at the endpoints.  To do that we add two fake data points, $\tg_0$ and
$\tg_{N+1}$ outside our interval and then we get:
\begin{alignat}{2}
\frac{1}{6}\left(\tg_{i-1}+4\tg_1+\tg_{i+1}\right) 
&= g_i  &\quad&i = 1,\ldots,N. 
\end{alignat}
And for the auxiliary condition we can either use a clamped spline which
gives the constraint:
\begin{eqnarray}
  \frac{1}{2}\left(-\tg_{0}+\tg_{2}\right) &=& \Delta g'_{1}
  \nonumber \\
  \frac{1}{2}\left(-\tg_{N-1}+\tg_{N+1}\right) &=& \Delta g'_{N}
\end{eqnarray}
or we can use a natural spline which gives the constraint:
\begin{eqnarray}
  \tg_{0}-2\tg_{1}+\tg_{2} &=& 0
  \nonumber \\
  \tg_{N-1}-2\tg_{N}+\tg_{N+1} &=& 0.
\end{eqnarray}
It is easiest to solve for $\tg_{0}$ and $\tg_{N+1}$ first and then put those
solutions into the equations we had before.  Then we obtain in the clamped
case, here on matrix form:
\begin{equation}
\frac{1}{6}
\begin{pmatrix}
  4      & 2      & 0       & \ldots & 0      \\
  1      & 4      & 1       & \ddots & \vdots \\
  0      & \ddots & \ddots  & \ddots & 0      \\
  \vdots & \ddots &    1    &   4    & 1      \\
  0      & \ldots &    0    &   2    & 4      \\
\end{pmatrix}
\begin{pmatrix}
  \tg_{1} \\
  \tg_{2}  \\
  \vdots   \\
  \tg_{N-1} \\
  \tg_{N} \\
\end{pmatrix}
= 
\begin{pmatrix}
  g_{1} + \frac{1}{3}\Delta g'_{1} \\
  g_{2}   \\
  \vdots  \\
  g_{N-1} \\
  g_{N} + \frac{1}{3}\Delta g'_{N} \\
\end{pmatrix}
\end{equation}
and in the natural case
\begin{equation}
\frac{1}{6}
\begin{pmatrix}
  6      & 0      & 0       & \ldots & 0      \\
  1      & 4      & 1       & \ddots & \vdots \\
  0      & \ddots & \ddots  & \ddots & 0      \\
  \vdots & \ddots &    1    &   4    & 1      \\
  0      & \ldots &    0    &   0    & 6      \\
\end{pmatrix}
\begin{pmatrix}
  \tg_{1} \\
  \tg_{2}  \\
  \vdots   \\
  \tg_{N-1} \\
  \tg_{N} \\
\end{pmatrix}
= 
\begin{pmatrix}
  g_{1}   \\
  g_{2}  \\
  \vdots   \\
  g_{N-1} \\
  g_{N}   \\
\end{pmatrix}.
\end{equation}
With the auxiliary points the function $g$ is written as:
\begin{equation}
g(x) \simeq \sum_{i = 0}^{N+1}\tg_{i}\phi_{3}
\left(\frac{x-x_i}{\Delta}\right),
\end{equation}
so in the end we must not forget to compute the $\tg$'s for the auxiliary
points also.  We note here also that if we wanted to use higher order splines
we would need to add more auxiliary points outside the original dataset.

\section{Solution of the linear system}
The linear equations in the previous sections can be solved very efficiently
since the matrices involved are tridiagonal.  Denoting the matrix by $A$ we
solve the system by using $LU$ decomposition writing $A = LU$.  We will
attempt a solution of the following form:
\begin{align}
&
\begin{pmatrix}
  a_{11}   & a_{12} & 0       & \ldots    & 0      \\
  a_{21}   & a_{22} & \ddots  & \ddots    & \vdots \\
  0        & \ddots & \ddots  &           & 0      \\
  \vdots   & \ddots &         &           & a_{n-1,n}      \\
  0        &        &    0    & a_{n,n-1} & a_{n,n}      \\
\end{pmatrix} 
\quad = 
\nonumber \\ 
&
\begin{pmatrix}
  l_{11}   & 0      & 0       & \ldots    & 0      \\
  l_{21}   & l_{22} & \ddots  & \ddots    & \vdots \\
  0        & \ddots & \ddots  &           &        \\
  \vdots   & \ddots &         &           & 0              \\
  0        &        &    0    & l_{n,n-1} & l_{n,n}      \\
\end{pmatrix}
\begin{pmatrix}
  1        & u_{12} & 0       & \ldots    & 0      \\
  0        & \ddots & \ddots  & \ddots    & \vdots \\
  \vdots   & \ddots &         &           &        \\
           &        &         &           & u_{n-1,n}  \\
  0        & \ldots &         &     0     & 1           \\
\end{pmatrix}.
\end{align}
Once we have determined $L$ and $U$ we can solve the equation $Ax = b$ by
first solving $Lz = b$ and then $Ux = z$ both of which are trivial to solve.
\begin{alignat}{2}
a_{11}    &= l_{11}    &&
\nonumber \\
a_{i,i-1} &= l_{i,i-1}  
&\qquad&\mbox{for each } i = 2,3,\ldots,N;
\nonumber \\
a_{i,i} &= l_{i,i-1}u_{i-1,i}+l_{i,i}  
&\qquad&\mbox{for each } i = 2,3,\ldots,N;
\nonumber \\
a_{i,i+1} &= l_{i,i}u_{i,i+1}
&\qquad&\mbox{for each } i = 2,3,\ldots,N.
\end{alignat}
The solution is found in the following manner:
\begin{enumerate}
\item Step 1: Set 
$$
l_{11} = a_{11}
\qquad\mbox{and}\qquad 
u_{12} = a_{12}/l_{11}
$$
\item Step 2: For $i = 2,\ldots,N-1$ set
$$
l_{i,i-1} = a_{i,i-1}, 
\qquad
l_{i,i} = a_{i,i}-l_{i,i-1}u_{i-1,i} 
\qquad\mbox{and}\qquad 
u_{i,i+1} = a_{i,i+1}/l_{i,i}
$$
\item Step 3: Set
$$
l_{n,n-1} = a_{n,n-1}
\qquad\mbox{and}\qquad 
l_{n,n} = a_{n,n}-l_{n,n-1}u_{n-1,n}.
$$
\end{enumerate}


\section{Integration of splines and endpoint corrections}
Once the linear coefficients $\tg_{i}$ are found the function can be
evaluated in terms of the B-splines and we can for instance use the
B-spline representation to obtain accurate integrals and transforms of
the original function.  Let us consider how that can be done for the
cubic B-spline case.  We consider the integral of $g$ with a kernel function
$K(x)$ that we assume that is defined on the whole real axis:
\begin{eqnarray}
\int_{x_1}^{x_N}f(x)g(x)dx &\simeq& 
\sum_{i = 0}^{N+1}\tg_{i}\int_{x_1}^{x_N}K(x)
\phi_{3}\left(\frac{x-x_i}{\Delta}\right)
\nonumber\\
&=& 
\sum_{i = 0}^{N+1}\tg_{i}
\Delta\int_{\frac{x_1-x_i}{\Delta}}^{\frac{x_N-x_i}{\Delta}}
K(\Delta~x+x_i)\phi_{3}(x)dx
\nonumber\\
&=& 
\sum_{i = 0}^{N+1}\tg_{i}
\Delta\int_{1-i}^{N-i}
K(\Delta~x+x_i)\phi_{3}(x)dx
\nonumber\\
&=& 
\sum_{i = 1}^{N}\tg_{i}K_{i}+
\Bigl[ 
\tg_{0}\kappa_{0}
+\tg_{N+1}\kappa_{N+1}
\Bigr]
\nonumber\\ &-&
\Bigl[
\tg_{1}\kappa_{1}
+\tg_{N}\kappa_{N}
+\tg_{2}\kappa_{2}
+\tg_{N-1}\kappa_{N-1}
\Bigr],
\end{eqnarray}
where we have defined the following quantities:
\begin{alignat}{2}
  &K_{i} = \Delta\int_{-\infty}^{\infty}K(\Delta~x+x_i)\phi_{3}(x)dx 
  & &
  \nonumber \\
  &\kappa_{0} = \Delta\int_{1}^{2}K(\Delta~x+x_0)\phi_{3}(x)dx
  &\qquad&
  \kappa_{N+1} = \Delta\int_{-2}^{-1}K(\Delta~x+x_{N+1})\phi_{3}(x)dx
  \nonumber \\
  &\kappa_{1} = \Delta\int_{-2}^{0}K(\Delta~x+x_1)\phi_{3}(x)dx
  &\qquad&
  \kappa_{N} = \Delta\int_{0}^{2}K(\Delta~x+x_{N})\phi_{3}(x)dx
  \nonumber \\
  &\kappa_{2} = \Delta\int_{-2}^{-1}K(\Delta~x+x_2)\phi_{3}(x)dx
  &\qquad&
  \kappa_{N-1} = \Delta\int_{1}^{2}K(\Delta~x+x_{N-1})\phi_{3}(x)dx
\end{alignat}

In the final expression we have split the result up into two parts;
one part is a simple sum that has just $N$ datapoints and where each
integral has the same limits, i.e. from -2 to +2, here written from
$-\infty$ to $\infty$ for simplicity.  The second part is the 
endpoint correction and is the contribution from the auxiliary
datapoints and a correction for the change in the integration limits
that we introduced into the first term.  We note that in case we are
using other B-splines such as the linear or higher order ones, the
first term would have the same form whereas the endpoint corrections
would be different.

\section{Hilbert Transforms}
Hilbert transforms come up quite frequently in physics whenever one is dealing
with response functions of some sort describing the reaction of physical
systems to external influence.  Typical examples of such response functions
are dielectric functions in electrodynamics and Green's function and self
energies one encounters in quantum many body physics.  Causality will require
such functions to be analytic as a function of frequency in one half plane and
their real and imaginary parts are related by the celebrated Kramers-Kronig
relations.  Quite often one only has knowledge of either the real part or the
imaginary part of these function but not both and thus one is forced to
calculate a Hilbert transform in order to obtain the other half.

Let us consider a real function $g(x)$.  We will assume that $g$ is of
essentially bounded support, i.e. it is zero for frequencies of
magnitude greater than some given value.  In addition we assume that
$g$ is non-negative as well.  This last assumption is not necessary
but is quite frequently true in practice.  The Hilbert transform of
$g$ is now defined as follows;
\begin{equation}
\hat{g}(y) = \principal_{-\infty}^{\infty}{g(x)dx \over y-x}
\end{equation}
Here P denotes Cauchy principal value.  We prefer to extend this
definition a little bit and define for a complex $y$:
\begin{equation}
\hat{g}(y) = \int_{-\infty}^{\infty}{g(x)dx \over y-x}
\end{equation}
and if $\imag y = 0$ we compute:
\begin{equation}
\hat{g}(y) = \int_{-\infty}^{\infty}{g(x)dx \over y+i0-x}
= \principal_{-\infty}^{\infty}{g(x)dx \over y-x}-i\pi g(y).
\end{equation}
Let us now compute this Hilbert transform using the cubic B-spline.
The integration kernel function is here:
\begin{equation}
K(y,x) = \frac{1}{y-x}
\end{equation}
and thus we get:
\begin{alignat}{2}
  &K_i(y) = \int_{-\infty}^{\infty}
  \frac{1}{\frac{y-x_{i}}{\Delta}-x}\phi_{3}(x)dx 
  & &
  \nonumber \\
  &\kappa_{0}(y)  = \int_{1}^{2}
  \frac{1}{\frac{y-x_{0}}{\Delta}-x}\phi_{3}(x)dx
  &\qquad&
  \kappa_{N+1}(y) = \int_{-2}^{-1}
  \frac{1}{\frac{y-x_{N+1}}{\Delta}-x}\phi_{3}(x)dx
  \nonumber \\
  &\kappa_{1}(y)  = \int_{-2}^{0}
  \frac{1}{\frac{y-x_{1}}{\Delta}-x}\phi_{3}(x)dx
  &\qquad&
  \kappa_{N}(y)   = \int_{0}^{2}
  \frac{1}{\frac{y-x_{N}}{\Delta}-x}\phi_{3}(x)dx
  \nonumber \\
  &\kappa_{2}(y)  = \int_{-2}^{-1}
  \frac{1}{\frac{y-x_{2}}{\Delta}-x}\phi_{3}(x)dx
  &\qquad&
  \kappa_{N-1}(y) = \int_{1}^{2}
  \frac{1}{\frac{y-x_{N-1}}{\Delta}-x}\phi_{3}(x)dx.
\end{alignat}
It is useful to introduce the following functions:
\begin{equation}
  H(z) = \int_{-\infty}^{\infty}\frac{\phi_{3}(x)}{z-x}dx
  \qquad
  H_{1}(z) = \int_{0}^{2}\frac{\phi_{3}(x)}{z-x}dx
  \qquad
  H_{2}(z) = \int_{1}^{2}\frac{\phi_{3}(x)}{z-x}dx
\end{equation}
and in terms of these we have
\begin{alignat}{2}
  &K_{i}(y) = H\left(\frac{y-x_1}{\Delta}\right)
  & &
  \nonumber \\
  &\kappa_{0}(y)  = H_{2}\left(\frac{y-x_0}{\Delta}\right)
  &\qquad&
  \kappa_{N+1}(y) = -H_{2}\left(\frac{x_{N+1}-y}{\Delta}\right)
  \nonumber \\
  &\kappa_{1}(y)  = -H_{1}\left(\frac{x_{1}-y}{\Delta}\right)
  &\qquad&
  \kappa_{N}(y)   = H_{1}\left(\frac{y-x_{N}}{\Delta}\right)
  \nonumber \\
  &\kappa_{2}(y)  = -H_{2}\left(\frac{x_{2}-y}{\Delta}\right)
  &\qquad&
  \kappa_{N-1}(y) = H_{2}\left(\frac{y-x_{N-1}}{\Delta}\right)
\end{alignat}

We now simply need to compute the Hilbert transforms for the $H$ functions
which is easily done in Mathematica.
\begin{eqnarray}
  H(y) &=& \frac{1}{6}
        \Bigl\{
        l_{3}(y+2)-4l_{3}(y+1) + 6l_{3}(y) - 4l_{3}(y-1)+l_{3}(y-2)     
        \Bigr\}
  \nonumber \\
  H_{1}(y) &=& \frac{1}{6}
        \Bigl\{
        l_{3}(y-2)-4l_{3}(y-1)
        +3l_{3}(y)-6l_{2}(y)
        \Bigr\}
        +\frac{11-3y^2}{9}+\frac{2}{3}\log(-y)  
 \nonumber \\
  H_{2}(y) &=& \frac{1}{6}
        \Bigr\{
        l_{3}(2-y)-l_{3}(y-1)+
        3l_{2}(y-1)-3l_{1}(y-1)
        \Bigr\}
  \nonumber \\ &&       
        +\frac{32-27y+6y^2}{36}+\frac{1}{6}\log(1-y).
\end{eqnarray}
Here we have introduced the function $l_{n}(x) \equiv x^{n}\log(-x)$
and separated the divergent term of the transform from the regular term.

For large value of $y$ it is useful to have asymptotic expansion in inverse
powers of $y$. For $H$ this expansion is quite simple since $H$ is an odd
function:
\begin{equation}
  H(y) = 
  \frac{1}{y}
  \left\{1+\frac{1}{y^2}
    \left\{\frac{1}{3}+\frac{1}{y^2}
      \left\{\frac{3}{10}+\frac{1}{y^2}
        \left\{\frac{17}{42}+\frac{1}{y^2}
          \left\{\frac{32}{45}+\frac{1}{y^2}\frac{15}{11}
          \right\}
        \right\}
      \right\}
    \right\}
  \right\}.
\end{equation}
For the other ones we write the expansion as:
\begin{equation}
  h_0+\frac{1}{y}
  \left\{h_1+\frac{1}{y}
    \left\{h_2+\frac{1}{y}
      \left\{h_3+\frac{1}{y}
        \Biggl\{h_4+\ldots
        \Biggr\}
      \right\}
    \right\}
  \right\},
\end{equation}
where the expansion coefficients are:
\begin{gather}
  H_1: \quad 
  \left\{
    0,\frac{1}{2},\frac{7}{30},\frac{1}{6},\frac{31}{210},\frac{3}{20},
    \frac{127}{756},\frac{17}{84},\frac{511}{1980},\frac{31}{90},
    \frac{2047}{4290},\frac{15}{22},\frac{8191}{8190},\ldots
  \right\}
  \nonumber \\
  H_2: \quad
  \left\{
    0,\frac{1}{24},\frac{1}{20},\frac{11}{180},\frac{8}{105},
    \frac{163}{1680},\frac{191}{1512},\frac{53}{315},\frac{227}{990},
    \frac{3797}{11880},\frac{3907}{8580},\frac{7957}{12012},
    \frac{4024}{4095},\ldots
  \right\}.
\end{gather}

If $y$ is complex we can use complex logarithms in our computation, however if
$y$ is real we should be a little bit more careful.  In that case the argument
in the logarithm should be replaced by its absolute value and we treat the
imaginary part separately.  We get an imaginary contribution only if $y$ is
somewhere inside the knots and the value is simply $-\pi\phi_{3}(y)$ as we can
easily see by looking at the equation above and considering the limit
$y+i\epsilon$ when $\epsilon$ goes to zero.

We should note that if $y$ coincides with one of the knots we should in
general simply omit the contribution from that knot since $x^{n}\log(x)$ goes
to zero as $x$ goes to zero for $n>0$.  However if $g$ is discontinuous at the
boundaries and $y$ happens to coincide with one of the knots that are used in
the endpoint corrections we will get logarithmic divergences in the transform
as is evident from the equations for the $H$ functions. In most cases we
encounter however the function $g$ goes continuously to zero at the boundary
and we have $g_1 = g_N = 0$.  It is easy to see from the spline equations that
the divergent terms in the endpoint corrections will cancel in this case.

We should remark here at the end that it is necessary to include the endpoint
corrections even though the integration limits are unbounded and the function
is zero at the endpoints of our numerical integration limits.  The reason for
this is that the spline is not valid outside $[x_1,x_N]$ and in particular will
not be equal to zero just outside this interval.


\section{Fourier Transforms}
Now we turn to Fourier transforms which as everyone knows are ubiquitous in
science and more useful and important than the Hilbert transforms.  In keeping
with notational convention we replace $x_{i}$ by $t_{i}$ and $y$ by $\omega$.
The kernel function is now
\begin{equation}
  K(\omega) = \exp(i\omega t)
\end{equation}
and it is easy to see that we get:
\begin{alignat}{2}
  &K_{i}(y) = \Delta\exp(i\omega t_{i})F(i\omega\Delta)
  & &
  \nonumber \\
  &\kappa_{0}(y)  = \Delta\exp(i\omega t_{0})F_{2}(i\omega\Delta)
  &\qquad&
  \kappa_{N+1}(y) = \Delta\exp(i\omega t_{N+1})\overline{F_{2}(i\omega\Delta)}
  \nonumber \\
  &\kappa_{1}(y)  = \Delta\exp(i\omega t_{1})\overline{F_{1}(i\omega\Delta)}
  &\qquad&
  \kappa_{N}(y)   = \Delta\exp(i\omega t_{N})F_{1}(i\omega\Delta)
  \nonumber \\
  &\kappa_{2}(y)  = \Delta\exp(i\omega t_{2})\overline{F_{2}(i\omega\Delta)}
  &\qquad&
  \kappa_{N-1}(y) = \Delta\exp(i\omega t_{N-1})F_{2}(i\omega\Delta)
\end{alignat}
where we have defined:
\begin{equation}
  F(z) = \int_{-\infty}^{\infty}e^{zt}\phi_{3}(x)dx
  \qquad
  F_{1}(z) = \int_{0}^{2}e^{zt}\phi_{3}(x)dx
  \qquad
  F_{2}(z) = \int_{1}^{2}e^{zt}\phi_{3}(x)dx
\end{equation}
These integrals are easy to do and we obtain:
\begin{eqnarray}
  F(z) &=& \left(\frac{\sinh(z/2)}{z/2}\right)^{4}
  \nonumber\\
  F_{1}(z) &=& \frac{1}{3z^{4}}\Bigl(3\exp(2z)-12\exp(z)+9+6z-2z^3\Bigr).
  \nonumber\\
  F_{2}(z) &=& \frac{\exp(z)}{6z^{4}}\Bigl(6\exp(z)-6-6z-3z^2-z^3\Bigr)
\end{eqnarray}
It is also useful here to have series expressions for these functions for
small values of $z$ and those are given by:
\begin{equation}
  F(z) = 1+y^2
  \left\{\frac{1}{6}+y^2
    \left\{\frac{1}{80}+y^2
      \left\{\frac{17}{30240}+y^2
        \left\{\frac{31}{1814400}+y^2
          \frac{1}{2661120}
        \right\}
      \right\}
    \right\}
  \right\}.
\end{equation}
The other ones we write on the form
\begin{equation}
  f_0+\frac{1}{y}
  \left\{f_1+\frac{1}{y}
    \left\{f_2+\frac{1}{y}
      \left\{f_3+\frac{1}{y}
        \Biggl\{f_4+\ldots
        \Biggr\}
      \right\}
    \right\}
  \right\},
\end{equation}
where the expansion coefficients are:
\begin{gather}
  F_{1}: 
  \left\{
    \frac{1}{2},\frac{7}{30},\frac{1}{12},\frac{31}{1260},
    \frac{1}{160},\frac{127}{90720},\frac{17}{60480},
    \frac{73}{1425600},\frac{31}{3628800},\frac{2047}{1556755200},
    \frac{1}{5322240}
  \right\}
  \nonumber\\
  F_{2}: 
  \left\{
    \frac{1}{24},\frac{1}{20},\frac{11}{360},\frac{4}{315},
    \frac{163}{40320},\frac{191}{181440},\frac{53}{226800},
  \right.
  \nonumber \\
  \left.        
    \frac{227}{4989600},\frac{3797}{479001600},
    \frac{3907}{3113510400},\frac{7957}{43589145600}
  \right\}.
\end{gather}
Let us now look more carefully at the first term in our expression for the
Fourier transform.  This is the term:
\begin{equation}
  I_{0} = \sum_{i=1}^{N}\tg_{i}\Delta\exp(i\omega t_{i})F(i\omega\Delta)
  = F(i\omega\Delta)\exp(i\omega x_{0})\Delta\sum_{n=0}^{N-1}\tg_{n+1}\exp(i\omega n\Delta).
\end{equation}
If we now take $\omega$ of the form,
\begin{equation}
  \omega_k = \frac{2\pi k}{N\Delta}
\end{equation}
then we get,
\begin{eqnarray}
  I_{0} &=& F(i2\pi k/N)\exp(i\omega_{k} x_{0})
  \Delta\sum_{n=0}^{N-1}\tg_{n+1}\exp(i2\pi kn/N) 
  \nonumber \\ &=& 
  F(i2\pi k/N)\exp(i\omega_{k} x_{0})
  \Delta\fft[\tg](k).
\end{eqnarray}
Here we have used the notation $\fft[\tg](k)$ to denote the Fast Fourier
Transform of the vector $(\tg_1,\ldots,\tg_N)$.  We should note here that even
though $\omega$ is taken to be of the particular form shown above there is no
limit on the value of the integer $k$ and no aliasing is present even though
we take $k$ larger than $N/2$.  This is due to the presence of the attenuation
factor $F(i2\pi k/N)$ which ensures that this expression returns the "exact"
value of the Fourier integral.

It is clear that if we are computing the Fourier transforms of non-periodic
functions the inverse Fourier transform can be treated in the same fashion.
For periodic functions, like Green's functions in complex time, the inverse
Fourier transform is a simple series and thus we do not need any integration
at all.

\section{Fourier transforms for Green's functions on Matsubara axis}
To take a concrete example let us now look at Fourier transforms for Green's
functions on the Matsubara axis.  Here we let $g$ be a Green's function which
can be either Fermionic or Bosonic defined on the Matsubara Frequencies
$\omega_{n} = n\pi T$.  If $g$ is Fermionic it is nonzero only on the odd
frequencies and on the even frequencies if it is Bosonic.  I also assume that
$\overline{g(i\omega_{n})} = g(-i\omega_{n})$ but we should note that this
relation does not necessarily hold for the off-diagonal elements of
multiorbital, matrix Green's function.  Now let
\begin{eqnarray}
x = {\tau \over \beta} &=& {k \over N/2}, \qquad k = 0,\ldots,N/2-1
\nonumber \\
 &=& {k \over N/2}-2, \qquad k = N/2,\ldots,N-1
\end{eqnarray}
Then we have
\begin{eqnarray}
g(\tau) &=&
T\sum_{n}e^{-i\omega_{n}\tau}g(i\omega_{n})
\nonumber \\
&=& T\sum_{n=1}^{\infty}e^{-i\pi nx}g(n)
+ T\sum_{-\infty}^{-1}e^{-i\pi nx}g(n)
+Tg(0)
\nonumber \\
&=& T\sum_{n=1}^{\infty}e^{-i\pi nx}g(n)
+ T\sum_{1}^{\infty}e^{-i\pi nx}g(-n)
+Tg(0)
\nonumber \\
&=& 2T \real \sum_{n=0}^{\infty}e^{-i\pi nx}g(n)
-Tg(0)
\nonumber \\
&\simeq& 2T \real\sum_{n=0}^{\infty}e^{-i2\pi nk/N}g(n)
-Tg(0)
\nonumber \\
&=& 2T \real\,\ifft[g](k) - Tg(0) 
\end{eqnarray}
Thus we see that we can use the inverse Fast Fourier Transform (iFFT) to get
the time dependent functions at discrete points, $k$.  

To go the other way we need to take into account the discontinuity of the
Green's Function at $\tau = m\times\beta$.  Here we spline the function on the
interval $[0^{+},\beta^{-}]$ and then the function can be obtained on the
interval $[-\beta^{+},0_{-}]$ using the (anti)periodicity property.  We should
note that the value of $g(\beta^{-})$ is not included in the data vector and
should be found from the known value of the discontinuity of $g$.

Define:
\begin{equation}
  \gamma_{i} = g(\tau_{i-1}) \quad i = 1,\ldots,N/2 \qquad\mbox{and}\qquad
  \gamma_{N/2+1} = g(\beta^{-}).
\end{equation}
Then $\left\{\tgam\right\}_{i=0}^{N/2+2}$ is the set of spline coefficients
for the function $g(\tau)$ in the interval $[0^{+},\beta^{-}]$ and
$\left\{\eta\tgam\right\}_{i=0}^{N/2+2}$ is the corresponding set for the
interval $[-\beta^{+},0^{-}]$ where $\eta = e^{i\omega_{n}\beta}$ is +1 for
bosons and -1 for fermions.  Define now also $\tg_{i} = \tgam_{i+1}$ for $i =
0,\ldots,N/2-1$ and $\tg_{i} = \eta\tgam_{i+1-N/2}$ for $i = N/2,\ldots,N-1$.
Then we get (using $\theta = \omega_{n}\Delta$),
\begin{eqnarray}
  g(i\omega_{n}) &=& \frac{1}{2}\int_{-\beta}^{\beta}
  e^{i\omega_{n}\tau}g(\tau)d\tau \quad = \quad 
  \frac{1}{2}\Delta F(i\theta)\sum_{k=0}^{N-1}e^{i\omega_{n}t_{k}}\tg_{k}
  \nonumber \\ &+& 
  \Delta\left\{
          e^{-i\theta}          F_{2}(i\theta) \tgam_{0}
    -                 \overline{F_{1}(i\theta)}\tgam_{1}
    -     e^{ i\theta}\overline{F_{2}(i\theta)}\tgam_{2}
  \right\}
  \nonumber \\ &+&
  \Delta\eta\left\{
          e^{ i\theta}\overline{F_{2}(i\theta)}\tgam_{N/2+2}
    +                 \overline{F_{1}(i\theta)}\tgam_{N/2+1}
    -     e^{-i\theta}          F_{2}(i\theta) \tgam_{N/2}
  \right\}
  \nonumber \\ &=& 
  \frac{1}{2}\Delta F(i\theta)\fft[\tg](n)
  +\Delta e^{-i\theta}F_{2}(i\theta)
  \Bigl[\tgam_{0}-\eta\tgam_{N/2}\Bigr] + 
  \nonumber \\ &+& 
  \Delta \overline{e^{i\theta}F_{2}}(i\theta)
  \Bigl[\eta\tgam_{N/2+2}-\tgam_{2}\Bigr] +
  \Delta \overline{F_{1}}(i\theta)
  \Bigl[\eta\tgam_{N/2+1}-\tgam_{1}\Bigr] 
\end{eqnarray}
The first sum in this expression is simply the Fast Fourier Transform sum.
Note that here we have used the fact that the endpoint corrections from the
lower interval are the same as the ones from the upper interval and note also
the plus sign on the fifth term in the endpoint correction.  The reason for
this change in sign is that here we have not included this point in the main
sum and therefore we have to add it's contribution.  Note also that since
$F_{2}$ always comes with $e^{-i\theta}$ tagging along we can combine them
beforehand and get a simpler expression for $F_{2}$ than we already have.

In order to check the validity of the formulas we should check to make sure
that they produce the correct asymptotic behavior.  Up to second order we
have,
\begin{equation}
  \Delta F_{1}(i\theta) \sim \frac{-2}{3i\omega}
  \qquad\mbox{and}\qquad
  \Delta e^{-i\theta}F_{2}(i\theta) \sim 
  \frac{-1}{6i\omega}+\frac{1}{2\Delta}\frac{1}{\omega^{2}}
\end{equation}
and thus,
\begin{eqnarray}
  g(i\omega) &\sim& \frac{1}{i\omega}
  \left(
    \frac{-\tgam_{0}}{6}+\frac{\eta\tgam_{N/2}}{6}+
    \frac{\eta\tgam_{N/2+2}}{6}+\frac{-\tgam_{2}}{6}+
    \frac{2\eta\tgam_{N/2+1}}{3}+\frac{-2\eta\tgam_{1}}{3}
  \right)
  \nonumber \\ && +
  \frac{1}{\omega^{2}}
  \left(
    \frac{\tgam_{0}-\eta\tgam_{N/2}}{2\Delta}+
    \frac{\eta\tgam_{N/2+2}-\tgam_{2}}{2\Delta}
  \right)
  \nonumber \\ &=&
  \frac{1}{i\omega}
  \left(
    \eta\frac{\tgam_{N/2}+4\tgam_{N/2+1}+\tgam_{N/2+2}}{6}-
    \frac{\tgam_{0}+4\tgam_{1}+\tgam_{2}}{6}
  \right)
  \nonumber\\ && +
  \frac{1}{\omega^{2}}
  \left(
    \eta\frac{\tgam_{N/2+2}-\tgam_{N/2}}{2\Delta}-
    \frac{\tgam_{2}-\tgam_{0}}{2\Delta}
  \right)
  \nonumber \\ &=&
  \frac{g(\tau = 0^{-})-g(\tau = 0^{+})}{i\omega}+
  \frac{g'(\tau = 0^{-})-g'(\tau = 0^{+})}{\omega^{2}}
\end{eqnarray}
Hence we see that the endpoint corrections reproduce correctly the first two
terms in the asymptotic expansion for the Green's function on the Matsubara
axis. 

\section{Other interpolation methods for Fourier Transforms}
We should here at the end mention that in Numerical Recipes \cite{Press:NR}
the computation of Fourier integrals using splines is somewhat different from
the presentation here.  The linear case, the formulas of which we have not
derived here, is exactly the same but in the cubic case some other basis
functions have been used for the spline.  The basis that is used there has the
benefit of not requiring the solution of a linear system since the basis
function is zero at all the knots except for the center one, see figure
\ref{fig:NRBasis}, but on the other hand it has the drawback of having
discontinuous derivatives at some or all of the knots.  This can be deduced
from the asymptotic behavior of the attenuation coefficient; the one derived
here falls off as $1/\omega^{4}$ indicating that the first two derivatives are
continuous, whereas the one derived in Numerical Recipes falls off as
$1/\omega^{2}$.  It does seem like a shame going through so much trouble in
using cubic interpolation for each subinterval and not having continuous
second derivatives over the whole integration region.

\begin{figure}[htb]
\centering
\begin{tabular}{c@{\;}c}
\mbox{\epsfig{file=Figures-App/BsNR1.ps,angle=-90,width=\halfpage}} &
\mbox{\epsfig{file=Figures-App/BsNR2.ps,angle=-90,width=\halfpage}} \\
\end{tabular}
\caption{Examples of basis functions of the sort used by Numerical
  Recipes. The figure on the right uses the same interpolation
  parameters as figure \ref{fig:B-splines}.  The basis function shown
  here can be written in terms of the cubic B-splines as
  $\left\{\right.$-3$\phi_{3}(x+3)$+12$\phi_{3}(x+2)$-45$\phi_{3}(x+1)$+168$\phi_{3}(x)$-45$\phi_{3}(x-1)$+12$\phi_{3}(x-2)$-3$\left.\phi_{3}(x-3)\right\}/97$.}
\label{fig:NRBasis}
\end{figure}   







%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
