Sinc Interpolation based on FFT

[stevengj]

@RainerHeintzmann

Maybe I can do a simulation of 100s of noisy signals which are filtered, to demonstrate how you probably get get a better estimate by the sqrt(2) scaling at the Nyquist sampling limit than by the mean (or sum as you mentioned).

You don’t need to do simulations, you can do this analytically. It simplifies a lot because (a) the DFT is linear and (b) the frequency components are orthogonal (when integrated over a common period), and the upshot is that you can analyze just the Nyquist term by itself:

In particular, take a function f(x)=a_{+}e^{i\pi x}+a_{-}e^{-i\pi x}, sampled at the Nyquist frequency f(n)=(a_{+}+a_{-})(-1)^{n} so we can only measure 2a=a_{+}+a_{-}. We reconstruct/interpolate it as a signal \tilde{f}(x)=2a\left[c_{+}e^{i\pi x}+c_{-}e^{-i\pi x}\right]. Question: what coefficients c_{\pm} minimize the expected mean-square error E[\int|f(x)-\tilde{f}(x)|^{2}dx] if a_{\pm} are i.i.d. random numbers with zero mean and some distribution (e.g. Gaussian)?

By explicit computation:

E\left[\frac{1}{2}\int_{0}^{2}|f(x)-\tilde{f}(x)|^{2}dx\right] \\ = \frac{1}{2} E\left[\int_{0}^{2}\left\{ |a_{+}-2ac_{+}|^{2}+|a_{-}-2ac_{-}|^{2}+\cancel{\#e^{2\pi ix}+\overline{\#}e^{-2\pi ix}}\right\} dx\right] \\ =E\left[|a_{+}-2ac_{+}|^{2}+|a_{-}-2ac_{-}|^{2}\right] \\ =E\left[|(1-c_{+})a_{+}|^{2}+|c_{+}a_{-}|^{2}+|(1-c_{-})a_{-}|^{2}+|c_{-}a_{+}|^{2}\right] \\ =E[|a_{\pm}|^{2}]\left(|1-c_{+}|^{2}+|c_{+}|^{2}+|1-c_{-}|^{2}+|c_{-}|^{2}\right) \\ =E[|a_{\pm}|^{2}]\left(1+2\left|\frac{1}{2}-c_{+}\right|^{2}+2\left|\frac{1}{2}-c_{-}\right|^{2}\right)

(using the facts that E[a_{+}\overline{a_{-}}]=0 and E[|a_{+}|^{2}]=E[|a_{-}|^{2}]), which is minimized for \boxed{c_{\pm}=\frac{1}{2}}.

That is, the optimal interpolant is \boxed{\tilde{f}(x)=a\left[e^{i\pi x}+e^{-i\pi x}\right]}, which corresponds to splitting the Nyquist amplitude 2a equally between the positive- and negative-frequency terms. This also coincides with the minimal mean-square slope interpolant from above, and has the nice property that it interpolates real signals f (a_- = \overline{a_+} \implies a purely real) with real interpolants \tilde{f}. (The same analysis also works, with the same optimal interpolant \tilde{f}, if we average over purely real signals f with random a_+=\overline{a_-}, since in that case we still have E[a_{+}\overline{a_{-}}]=E[a_+^2]=0 if the real and imaginary parts of a_+ are i.i.d. with zero mean.)

No \sqrt{2}.