Introduction
Radial basis function approximation has the form
$$f(x) \approx \sum_{j=1}^N c_j \phi(\| x-x_j\|)$$
for some function $\phi(r)$ such as a Gaussian.
Solving for the coefficients $c_j$ is a linear problem and might be done either as a least-squares fit or as an interpolation problem.
As mentioned on the Community posting, http://community.wolfram.com/groups/-/m/t/1255676, there is a basis for RBF interpolation over a regularly spaced grid that has nice properties: it has the "cardinal basis" property, is a numerically accurate approximation of Gaussian RBF, and, for the OP's purpose, has a known Fourier transform.
I am not an expert on RBF approximation, but it is known that RBF has some issues: "error saturation," in which the interpolation error does not vanish as the grid spacing approaches zero; and "boundary layer" at the end points of a truncated, finite interval. Some appreciation of error should be considered in the choice of a method and its parameters.
We will give code for a straightforward implementation of RBF interpolation over a regularly spaced grid over a truncated, finite interval. We next apply it to the OP's Airy function example. To give a nice interface to the Fourier transform of an RBF interpolation, we give a more sophisticated implementation of the interpolation, as well as its transform, and apply them to the OP's example. Finally, we illustrate the issues with error in RBF interpolation with regard to the OP's example. It will be seen that the RBF approximation is very good overall and worst in the neighborhood of the singularity of the HeavisideTheta factor at x == 0.
A cardinal basis for Gaussian radial basis function interpolation
A "cardinal basis" for a grid $x_j$ is a set of functions $C_j(x)$ such that
$$C_j(x_j)=1,\quad\text{and}\quad C_j(x_k)=0\ \text{for all $k\ne j$.}$$
The Lagrange interpolation basis is a cardinal basis for a finite grid. The particular advantage of such a basis is that the interpolation formula does not require solving an equation:
$$f(x) \approx \sum_j f(x_j)\,C_j(x)\,.\tag{1}$$
Maz'ya & Schmidt, Approximate Approximations, (2007)
and
Boyd & Wang, "An analytic approximation to the cardinal functions of Gaussian radial basis functions on an infinite lattice", (2009)
independently derived the basis for an infinite grid with spacing $h = x_j - x_{j-1}$
and shape parameter $\alpha$:
$$C_j(x)=C(x-x_j),\quad\text{where}\quad
C(x)= {\alpha ^2\over\pi}\,
{\sin \left({\pi x}/{h}\right) \over {\sinh}\left({\alpha^2 x}/{h}\right)}$$
Since for this grid, the formula in (1) is an infinite sum, we will have to truncate it and sum over a finite interval.
A basic implementation for this idea is the following (note that the sum (1) is effected with Dot):
(* cBasis[x_, α_, h_:1] := α^2/Pi Sin[Pi x/h]/Sinh[α^2 x/h]; (* risks divide-by-zero *) *)
cBasis[x_, α_, h_: 1] := Sinc[Pi x/h]/Sinc[I α^2 x/h]; (* uses complex numbers *)
rbfInterpolation[f_, {x_, x1_, x2_}, α_, h_] := (* interpolates f over x1 <= x <= x2 *)
With[{fvals = f /. x -> # & /@ Range[x1, x2, h]},
Function[x, fvals.cBasis[x - Range[x1, x2, h], α, h]]
];
The main problem is that cBasis[0, α, h] needs to evaluate to 1, and the first, commented code above yields Indeterminate (from 0/0). The second avoids dividing by zero but does its computations with complex numbers (with a zero imaginary part). One could define cBasis[x_ /; x == 0, α_, h_: 1] := 1 explicitly and make it Listable, but this is slower. In any case, I wanted to start with simpler code. A slightly improved (especially if x == 0 is rare), but somewhat advanced, version will be given below with a more convenient, if more complicated, rbfInterpolation. These improvements will make it easier to handle the Fourier transform of an RBF interpolation.
Application to OP's Airy function
ff[x_] := AiryAi[x + AiryAiZero[3]] UnitStep[x]; (* HeavsideTheta is undefined at 0 *)
testIF = rbfInterpolation[ff[x], {x, 0, 12}, 1/2, 1/8];
Here's a plot showing the OP's function (blue) covered by the interpolation function (gold), together with the logarithm base-10 of the error (green) below, which we plot beyond the interpolation range:
Plot[{20 ff[x], 20 testIF[x], ff[x] - testIF[x] // RealExponent},
{x, -2, 15},
Frame -> True, FrameTicks -> {
{Automatic, Charting`ScaledTicks[{20 # &, #/20 &}]},
{Automatic, Automatic}},
FrameLabel -> {{"Log error", HoldForm@f[x]}, {HoldForm[x], None}}]
An understanding of the error in RBF interpolation can help choosing α and h and help explain why the error decreases between the end points. See Error analysis below.
Improved interpolation code and Fourier transform
One concern is speed, only because this method already appears somewhat slow to me. When possible, the code should preserve packed arrays, and use arrays of real numbers for real values. Also, the Fourier transform of an RBF interpolation is simply
$$\hat{f}(t) \approx \sum_j f(x_j) \, \hat{C}_j(t)$$
where $\hat{C}_j(t) = e^{i x_j t} \hat{C}(t)$ and $\hat{C}$ is the Fourier transform of $C$, is given by
$$\hat{C}(t) = -\frac{h \, e^{i x_j t} \left(\tanh
\left(\frac{\pi (h t-\pi )}{2 \alpha^2}\right)-\tanh \left(\frac{\pi (h t+\pi
)}{2 \alpha ^2}\right)\right)}{2 \sqrt{2 \pi
}}\,.$$
The only data needed are the grid $x_j$ and the function values $f(x_j)$ on the grid. We can store these in a simple data structure, similar to InterpolatingFunction, instead of in a Function as above.
To fix the divide by zero problem, we can adapt @Carl Woll's answer, which has the property of preserving packed arrays, when the arguments do not lead to 0/0. Even when they do, the unpacking is faster than making cBasis have the attribute Listable.
ClearAll[cBasis]; (* cardinal basis *)
(* for symbolic analysis *)
cBasis /: Normal[cBasis[x_, α_, h_: 1]] := α^2/Pi Sin[Pi x/h]/Sinh[α^2 x/h];
division[a_, b_] := Quiet@Block[{Indeterminate = 1}, Divide[a, b]];
cBasis[x_?NumericQ | x_?(VectorQ[#, NumericQ] &), α_, h_: 1] :=
division[α^2/Pi Sin[Pi x/h], Sinh[α^2 x/h]];
ClearAll[rbfIF, rbfInterpolation]; (* RBF interpolation *)
rbfInterpolation[f_, {x_, x1_, x2_}, α_, h_] :=
Module[{xdata, fdata},
xdata = Range[x1, x2, h];
fdata = Developer`ToPackedArray[f /. x -> # & /@ xdata];
rbfIF[xdata, fdata, {α, h}]
];
rbfIF[xdata_, fdata_, {α_, h_}][x_] := fdata.cBasis[x - xdata, α, h];
rbfIF /: MakeBoxes[if_rbfIF, form_] := BoxForm`ArrangeSummaryBox[ (* optional formatting *)
rbfIF, if, None, {{MinMax@First@if}}, {}, form, "Interpretable" -> Automatic];
ClearAll[ft, rbfFT]; (* Fourier transform of the cardinal basis/RBF interpolation *)
Block[{k, x0, α, h, x},
ft[k_, x0_, α_, h_: 1] = (* pre-calculate basic Fourier transform (takes a few secs.) *)
Simplify[
Exp[I k x0] FourierTransform[Normal[cBasis[x, α, h]], x, k,
Assumptions -> α > 0 && h > 0]]
];
rbfFT[rbfIF[args___]] := rbfFT[args]; (* for the transform of a rbfIF, change the head *)
rbfFT[xdata_, fdata_, {α_, h_}][k_] := fdata.ft[k, xdata, α, h];
rbfFT /: MakeBoxes[ft_rbfFT, form_] := BoxForm`ArrangeSummaryBox[ (* optional formatting *)
rbfFT, ft, None, {{MinMax@First@ft}}, {}, form, "Interpretable" -> Automatic];
The expressions for cBasis and its Fourier transform are simple elementary expressions. For machine-precision input, they will be vectorized. Note that the symbolic form of cBasis is used to compute the expression for the Fourier transform. One might use different user-interfaces for the symbolic cBasis.
Fourier transform of the OP's Airy function
testIF = rbfInterpolation[ff[x], {x, 0, 15}, 1/2, 1/8]; (* re-evaluate with new defs *)
testFT = rbfFT@testIF; (* Fourier transform *)
For comparison here is a numerical integration for the transform (Integrate/FourierTransform failed):
ClearAll[nFT];
Module[{t},
li = NIntegrate`LevinIntegrandReduce[ff[x] Exp[I t x] /. UnitStep[_] :> 1, x];
With[{lirules = Sequence @@ (DeleteCases[li["Rules"], _["Variables", _]] /.
HoldPattern["DifferentialMatrices" -> dm_] :> "DifferentialMatrix" -> First[dm])},
nFT[t0_?NumericQ] := Block[{t = t0},
1/Sqrt[2 Pi] NIntegrate[
ff[x] Exp[I t x] /. {UnitStep[_] :> 1} // Evaluate,
{x, 0, Infinity},
Method -> {"LevinRule", lirules}]
]]];
The magnitude and argument of the transform of ff and the interpolation overlay each other:
Plot[{nFT[t], testFT[t]} // Abs // Evaluate, {t, -10, 10}] (* Compare Abs *)
Plot[{nFT[t], testFT[t]} // Arg // Evaluate, {t, -10, 10}, (* Compare Arg *)
MaxRecursion -> 3, Exclusions -> None]
Error analysis and the boundary layer phenomenon
There are at least three parameters affecting error, $\alpha$, $h$, and the "round-off plateau," which depends on working precision and is ordinarily around $\epsilon\,\|f\|_\infty$ or a little higher, where $\epsilon$ is the floating-point precision limit (e.g., $MachineEpsilon for machine precision). Connected with $\alpha$ and $h$, is the "boundary layer" phenomenon (see below), which affects the error near the end points of the interval.
Away from the boundary, the error seems to depend principally on $\alpha$ and $h$. According to Boyd, "Error saturation in Gaussian radial basis functions on a finite interval" (2010), the error tends to get worse if $\alpha$ gets too small or too large. Just by testing $\alpha = 0.4, 0.5, 0.6$ etc., I found $\alpha = 0.5$ to be good on the OP's example, although it varies slightly with $h$.
The boundary layer effect decreases as ~Exp[-α^2 x/h] according to
Boyd, op. cit. (2010). This can be seen for various $h$ in the dashed red envelopes in the figure below. The singular point at $x=0$ impose by the UnitStep/HeavisideTheta factor means that $x = 0$ is a boundary point no matter how far we extend the interpolation. It seems little can be done about the error near $x=0$. On the other hand one can see that at $x \approx 9$, the end point of the OP's interpolation example, the error is quite good for $h = 1/32, 1/64$, because the interpolation range was extended to $x = 12$, far enough out that the boundary layer effect is negligible compared to round-off. (The round-off plateau is around machine precision until the function becomes small and $x$ far enough away from the boundary-layer effect of the large peak near $x=4.5$; see the lower grid line at $\log_{10}\epsilon\approx-15.95$ in the figure below.)
testIF8 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./8];
testIF16 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./16];
testIF32 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./32];
testIF64 = rbfInterpolation[ff[x], {x, 0, 12}, 0.5, 1./64];
errplot = Plot[ (* errors of the approximations *)
ff[x] - {testIF8[x], testIF16[x], testIF32[x], testIF64[x]} // RealExponent // Evaluate,
{x, -2, 14},
Frame -> True, FrameLabel -> {HoldForm[x], "Log error"},
PlotStyle -> Opacity[0.7],
PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}],
PlotPoints -> 100]
boundarylayer = Show[ (* boundary layer envelopes *)
Plot[0.001 Range[5,2,-1] * (* constant factor found by eye *)
Exp[-0.5^2 (x)/(1./{8, 16, 32, 64})] // RealExponent // Evaluate,
{x, -2, 14},
PlotStyle -> Directive[AbsoluteThickness[1], Dashing[{0.02, 0.02}], Red]],
Plot[3*^-7 Exp[-0.5^2 (12 - x)/(1./{8, 16, 32, 64})] // RealExponent // Evaluate,
{x, -2, 14},
PlotStyle -> Directive[AbsoluteThickness[1], Dashing[{0.02, 0.02}], Red]]
];
Show[errplot, boundarylayer,
GridLines -> {None, {-6.6, -9.5, RealExponent[$MachineEpsilon]}},
Frame -> True,
PlotLabel -> Row[{"Boundary layer effects of ", HoldForm[f[x]], " at ",
HoldForm[x == {0, 12}], " (", HoldForm[α == 0.5], ", ",
HoldForm[h == {1/8, 1/16, 1/32, 1/64}], ")"}]]
Here is another view of the boundary layer effect at $x=9$ and what the inclusion of the basis functions beyond the boundary add to the interpolation. The basis functions, which in the interpolation are multiplied by the function values, contribute a relatively large amount to the interpolation near $x=9$. We show it for $h = 1/8$ instead of one of the more accurate grid spacings, because it is easier to see what is going on. One can compare the decreasing contribution below with the decreasing error above (see the top two grid lines).
With[{x1 = 9, x2 = 12, α0 = 0.5, h0 = 1./8},
Plot[
Table[cBasis[x - x0, α0, h0] // RealExponent, {x0, x1, x2, h0}] //Evaluate,
{x, x1 - 2, x2 + 2},
PlotStyle -> ReplacePart[
Table[Directive[AbsoluteThickness[1.], Opacity[0.5]], {x0, x1, x2, h0}],
{-1 -> Directive[Black, Opacity[1]]}],
GridLines -> {{x1}, Automatic}, GridLinesStyle -> Lighter@Blue,
Frame -> True,
FrameLabel -> {HoldForm[x], HoldForm[Log10@HoldForm[C][x - x0]]},
PlotLabel ->
Row[{"Boundary layer effects at ", HoldForm[x == x1], " (",
HoldForm[α == α0], ", ", HoldForm[h == h0], ")"}]]
]
The following show the magnitude of the error in the Fourier transform of ff and the approximation as well as the magnitude of the phase error, for the four approximations.
Plot[nFT[t] -
Through[(rbfFT @@@ {testIF8, testIF16, testIF32, testIF64})[t]] //
RealExponent // Evaluate,
{t, -10, 10},
MaxRecursion -> 3, Exclusions -> None,
PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}],
PlotLabel -> "Log error"]
Plot[Arg@nFT[t] -
Arg@Through[(rbfFT @@@ {testIF8, testIF16, testIF32, testIF64})[t]] //
RealExponent // Evaluate,
{t, -10, 10},
MaxRecursion -> 3, Exclusions -> None,
PlotLegends -> Thread[HoldForm[h] == 1/{8, 16, 32, 64}],
PlotLabel -> "Log error in Arg"]
FitandLinearModelFit,NonlinearModelFit, etc. would do the job. $\endgroup$