6
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • $\begingroup$ I was wondering if you have a citation for your footnote (1) in the other post. $\endgroup$
    – Victor Liu
    Aug 12, 2014 at 18:38
  • $\begingroup$ @VictorLiu: Not really, but it's something you can test easily. Basically, if you sample a function and FFT the samples and then compare the results to an analytic expression for the Fourier series coefficients (assuming an analytic expression exists), then the coefficients which are large will be very accurate (from a percent error perspective) compared to the analytic result, whereas the coefficients which are very small will be very bad, and the smaller they are, the worse they get. I can add an example which shows the effect, if you think it'd help the question. $\endgroup$
    – DumpsterDoofus
    Aug 12, 2014 at 20:05
  • $\begingroup$ I was asking because I'm encountering similar behavior for multi-dimensional FFTs. In my case, I'm not sure how much of a problem high relative error is if the absolute error is small. $\endgroup$
    – Victor Liu
    Aug 12, 2014 at 23:55
  • $\begingroup$ @VictorLiu: "I'm not sure how much of a problem high relative error is if the absolute error is small" It depends on the usage. Small absolute error guarantees excellent signal reconstruction by Parseval's theorem, and that's normally what most people are after, so it's perfectly good. However, if the purpose of the computation is to get small-modulus Fourier coefficients, then it's pretty bad. Part of the reason I went down this path in the first place was because someone had a question about band structure and I needed to compute high-frequency Fourier coefficients of a smooth surface, and $\endgroup$
    – DumpsterDoofus
    Aug 13, 2014 at 0:00
  • $\begingroup$ ...ran into exactly this problem. :) $\endgroup$
    – DumpsterDoofus
    Aug 13, 2014 at 0:01

1 Answer 1

Reset to default
3
$\begingroup$

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
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks, that makes sense; I guess I was just assuming that the left and right line segments to the center point were all that was needed, and zeroed the outer two segments, which (as you pointed out) is wrong. Do you have any idea as to how the $\varphi$ functions are calculated? $\endgroup$
    – DumpsterDoofus
    Aug 20, 2014 at 20:04
  • $\begingroup$ @DumpsterDoofus No idea; they don't say which segments of which Lagrange polynomials they use. For trapezoidal rule, $\varphi_0(s)$ is $-\psi(s)$ for $s<0$ and $0$ for $s>0$ because the goal is to cancel out the kernel's contribution from outside the domain of integration. I tried a few choices for cubic case but they didn't work. $\endgroup$
    – Kirill
    Aug 20, 2014 at 21:44
  • $\begingroup$ @DumpsterDoofus I edited the answer. $\endgroup$
    – Kirill
    Aug 20, 2014 at 22:23
  • $\begingroup$ @DumpsterDoofus They also don't say explicitly how they choose to enforce the condition $W(0)=1$. $\endgroup$
    – Kirill
    Aug 20, 2014 at 22:28
  • $\begingroup$ @DumpsterDoofus Page 699 of my edition: "To our knowledge, the cubic-order formulas derived here have not previously appeared in the literature." $\endgroup$
    – Kirill
    Aug 20, 2014 at 23:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.