THE PROGRAM

The algorithm suggested by the results of the previous section is to compute from the time Fourier transformed original data $U({\omega},x)$, for each frequency $\omega$, a prediction filter along the x axis which is then used at the frequency $n\times \omega$ to solve for the missing data. I have implemented such an algorithm in the case where n is equal to two.

Since the operations are performed in the $({\omega},x)$ domain, the data u(t,x) as well as a zero padded version u_pad(t,x) of the data, with twice as many time samples as u(t,x), are both Fourier transformed over the time axis. The one-step forward prediction filters $(P_{1}({\omega}),{\cdots},P_{L}({\omega}))$ are then computed for all the frequencies from zero to half the Nyquist by solving using Levinson's recursion the following linear system

$${\bf R}({\omega})\pmatrix{1\cr -P_{1}({\omega})\cr \vdots\cr -P_{L}({\omega})\cr}= \pmatrix{e\cr 0\cr \vdots\cr 0\cr}$$ (15)

where ${\bf R}({\omega})$ is the autocorrelation matrix of a frequency slice from $U_{pad}({\omega},x)$ and e is the minimum prediction error energy.

According to equation (12), for n=2, $(P_{1}({\omega}),{\cdots},P_{L}({\omega}))$ predicts one step ahead the interpolated data ${U'}(2{\omega},x)$. This relation may be thought of as a convolutional operation

$$P({\omega}){\ast}{U'}(2{\omega},x)=0$$ (16)

where $P({\omega})=(P_{1}{\omega},{\cdots},P_{L}({\omega}))$ and ${U'}(2{\omega},x)$ contains the original data interlaced with missing values. $P({\omega})$ is then separated into an even part and an odd part, respectively $P_{e}({\omega})$ and $P_{o}({\omega})$, and ${U'}(2{\omega},x)$ into missing and known data, respectively $M(2{\omega},x)$ and $K(2{\omega},x)$. $M(2{\omega},x)$ is then obtained from

$$M(2{\omega},x)=-P_{o}^{-1}({\omega}){\ast}(P_{e}({\omega}){\ast}K(2{\omega},x)).$$ (17)

The same computations are performed with the one-step backward prediction filters $(P_{1}^{\ast}({\omega}),{\cdots},P_{L}^{\ast}({\omega}))$ giving another set of values for the missing data. The two values found for the missing data in each case are then averaged to give the final result. The interpolated data ${U'}({\omega},x)$ is then time inverse fourier transformed to u'(t,x).

It is worth noticing that the backward prediction filters are easily computed from the forward prediction filters and that the zero padding in time of u(t,x) is necessary to avoid any type of filter interpolation in the $\omega$ direction between even frequencies $({\omega}_{k}=2k{\pi}/N$, with k even).