I asked a similar question on MathSE but with more added fluff, but didn't really get any straight answers, so I figured I'd ask here. Computing Fourier coefficients of a function using the FFT is notoriously inaccurate whenever the magnitude of the coefficients is small, and 13.9 of Numerical Recipes in C makes a partial fix which uses spline basis functions, allowing one to still use the FFT, but drastically reducing the number of samples required.The author uses trapezoidial and cubic corrections, and I was mainly interested in the cubic correction. I'm interested in the derivation because I want to see how it works for multidimensional corrections (like 2D functions).  The derivation of the $W(\theta)$ and $\alpha_j$ functions is not explicitly shown (since it's a code/recipe book), but the author states For linear interpolation (s) is piecewise linear, rises from 0 to 1  for s in (1, 0), and falls back to 0 for s in (0, 1). For higher-order interpolation, (s) is  made up piecewise of segments of Lagrange interpolation polynomials. It has discontinuous  derivatives at integer values of s, where the pieces join, because the set of points used in  the interpolation changes discretely.To that end, I tried to reproduce the cubic corrections by defining$$\psi(x)=\begin{cases}0&x<-1\\\text{Lag}_3\left(\begin{bmatrix}0&0&1&0\\-2&-1&0&1\end{bmatrix},x\right)&-1\leq x\leq0\\\text{Lag}_3\left(\begin{bmatrix}0&1&0&0\\-1&0&1&2\end{bmatrix},x\right)&0\leq x\leq1\\0&1<x\end{cases}$$where $$\text{Lag}_3\left(\begin{bmatrix}y_1&y_2&y_3&y_4\\x_1&x_2&x_3&x_4\end{bmatrix},x\right)=\sum_{j=1}^4y_k\prod_{\begin{smallmatrix}k=1\\ k\neq j\end{smallmatrix}}^4\frac{x-x_k}{x_j-x_k}$$is the cubic Lagrange polynomial which passes the points specified by the $x_k,y_k$. However, this yields (in Mathematica)CubicLagrangeInterpolation[X_, Y_, x_] :=   Sum[Y[[j]]*Product[If[j == k, 1, (x - X[[k]])/(X[[j]] - X[[k]])],            {k, 1, 4}], {j, 1, 4}]; \[Psi][x_] :=   Piecewise[{{0,      x < -1}, {CubicLagrangeInterpolation[{-2, -1, 0, 1}, {0, 0, 1,        0}, x], -1 <= x <= 0},          {CubicLagrangeInterpolation[{-1, 0, 1, 2}, {0, 1, 0, 0}, x],      0 <= x <= 1}, {0, x > 1}}]; FullSimplify[ Integrate[\[Psi][s]*Exp[I*\[Theta]*s], {s, -Infinity, Infinity}]]which yields$$W(\theta)=\int_{-\infty}^\infty\psi(s)e^{i\theta s}\,\mathrm{d}s=\frac{\theta ^2-2 \left(\theta ^2+3\right) \cos (\theta )-2 \theta  \sin (\theta )+6}{\theta ^4}$$which does not match the $W(\theta)$ provided in the article, although when plotted the two are somewhat similar in appearance. Likewise, I am confused how he obtains his $\varphi$ endpoint kernel functions. Does anyone know what the author meant when he said he uses Lagrange interpolation to obtain his result? Or am I simply using Lagrange interpolation incorrectly?Minor note: When I use linear Lagrange interpolation ($\text{Lag}_1$ functions) I get the correct result for $W(\theta)$ for the author's trapezoidial case, and $\psi(x)$ is just a unit triangle.
What spline functions are used in Section 13.9 of Numerical Recipes in C?
numerics;interpolation;fourier analysis;b spline
When I change your code to this:psi[x_] := Piecewise[  {{0, x < -2}   , {CubicLagrangeInterpolation[{-3, -2, -1, 0}, {0, 0, 0, 1}, x], -2 <= x <= -1}   , {CubicLagrangeInterpolation[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -1 <= x <= 0}   , {CubicLagrangeInterpolation[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1}   , {CubicLagrangeInterpolation[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}   , {0, x > 2}}]It gives me$$ \frac{8 \left(6+\theta ^2\right) \sin^4(\frac\theta2)}{3 \theta ^4}, $$which is equivalent to the NR formula$$ \frac{6+\theta^2}{3\theta^4}(3-4\cos\theta+\cos2\theta). $$I think what's going on here is that your kernel's domain is wrong; Lagrange polynomials might be zero on $\pm1,\pm2$, but they don't have to be zero in between, so long as $0$ is one of the interpolation nodes.Code for endpoint kernels.I can't quite understand why $\alpha_3$ is not zero, because the kernels only touch points at most a distance $3$ away, so shouldn't the point number $3$ be far away from the endpoints that there is no endpoint contribution? Note that $\alpha_3$ has different domain from $\alpha_{0,1,2},\psi$. They don't quite say which segments of which polynomials they use, so this is a bit confusing.Clear[lag]lag[X_, Y_, x_] /; Length[X] == Length[Y] :=  Sum[Y[[j]]*    Product[If[j == k, 1, (x - X[[k]])/(X[[j]] - X[[k]])], {k,       Length@X}], {j, Length@Y}];psi[x_] := Piecewise[   {{0, x < -2}    , {lag[{-3, -2, -1, 0}, {0, 0, 0, 1}, x], -2 <= x <= -1}    , {lag[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -1 <= x <= 0}    , {lag[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}    , {0, x > 2}}];Plot[psi[x], {x, -2, 2}]FullSimplify[ Integrate[psi[s] Exp[I*\[Theta]*s], {s, -Infinity, Infinity}]]alpha[0, x_] := Piecewise[   {{0, x < 0}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 0 <= x <= 1}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}    , {0, x > 2}}];alpha[1, x_] := Piecewise[   {{0, x < -2}    , {lag[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], -1 <= x <= 0}    , {lag[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}    , {0, x > 2}}];alpha[2, x_] := Piecewise[   {{0, x < -2}    , {lag[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -2 <= x <= -1}    , {lag[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -1 <= x <= 0}    , {lag[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}    , {0, x > 2}}];alpha[3, x_] := Piecewise[   {{0, x < -3}    , {lag[{-3, -2, -1, 0}, {0, 0, 0, 1}, x], -3 <= x <= -1}    , {lag[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -1 <= x <= 0}    , {lag[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1}    , {lag[{0, 1, 2, 3}, {1, 0, 0, 0}, x], 1 <= x <= 2}    , {0, x > 2}}];Plot[Evaluate@Table[alpha[k, s], {k, 0, 3}], {s, -2, 2}]Integrate[   Table[(alpha[k, s - k] - psi[s - k]) Exp[I*\[Theta]*s], {k, 0,      3}], {s, -Infinity, Infinity}] // ExpToTrig // FullSimplify{Series[Re[%] // ComplexExpand, {\[Theta], 0, 6}],    Series[Im[%] // ComplexExpand, {\[Theta], 0, 6}]} //   Transpose // TableForm