Tsai and Huang Frequency Domain Method (1D)

It is also easier if we begin our description with the one-dimensional case, as it eases the transition into the two-dimensional case of interest in image processing.  Where at all possible, we have maintained the identity of the variables used in the paper.
 

1)  Definitions

x         = real-space coordinate
w        = frequency-space coordinate
f(x)      = ideal image
p         = number of frames
   = ith sample in the kth frame (i will run from 0 to N-1)
N        = number of samples in each frame
T         = samling interval in space (i.e. f(x) is sampled every T units along x)
       = inital offset of kth frame.  Determined by motion registration.

F(w)    = Continuous Fourier Transform (CFT) of continuous signal f(x)
= CFT of shifted continuous signal 
     = nth component of the DFT of the kth frame of f(x) (i.e. nth component of DFT of )
      = sampling frequency (in spacial domain) = 1/T
L         = resolution improvement factor
J          = our complex number, not i, not j, because those are used extensively in this paper for other things.
 

2)  Method

We know from continuous Fourier Transform theory that there is a simple way to relate shifted functions to non-shifted functions:
             .
where F(w) is the continuous Fourier Transform of the ideal image f(x).  In images,however, we desire to look in discrete space, so we need a way to relate this to a disrete variable.  Given the discrete Fourier Transform,
             .
So we note here that n runs from 0 to N-1.  We then have p DFT vectors of length N, and since k{0,p-1} and n{0,N-1}.  We writeas the nth frequency from the kth frame - this now is essentially a 2D matrix.

Now we relate the DFT to the CFT with the "aliasing relationship":
             .
So we have now related each discrete component, that is the nth frequency component shifted by k, to the continuous representation shifted by k.  This becomes important.

Now is the part where we make the assumption that F(w) is bandlimited.  We do this by saying that F(w) = 0 for all w such that |w| > L * 

This enables us to write:
              .

where
              . which we note is a relation bewteen the DFT and the CFT
              . which comes from the combinations of exponentials

We notice that the matrix is something we already know.  We can calculate that with a DFT from the original data.  We notice also that  is a matrix that we can calculate simply knowing N, T, L, and all .

Note that  is a  matrix containing 2L points.  We have not yet determined the value of L.

Clearly now, what we wish to solve for is .  We can do this using a number of matrix techniques, which are described in the paper, but which we will not go into here in the 1D case.

Once we have solved this matrix equation, what do we do now?  We have values for this matrix , but what do these numbers mean?  By definition, they are the values of the continuous Fourier Transform at the point:
              .
So we have taken these 2L values and assigned w coordinates to them in the frequency domain as sample point of the frequency representation of the continuous ideal image.  Now again, the matrix equation above considers only the nth component of all p of our DFTs.  We have N-1 more matrix equations to solve, and thus 2L*(N-1) more points to add to our frequency representation of the continous ideal image.  In total, we will have 2LN points in our representation, whereas we started out with N (being the most points we could obtain from any one image.)  So we see that our resolution has increased by a factor 2L!  If we want a resolution increase of 2, we make L=1.  In general, we will use small-ish values of L, so perhaps the picture of the matrix equation above is misleading, but it was purposely drawn that way to illustrate the matrix multiplication relationships.

Of course, in order to recover the exact image exactly, we would need enough points to be above the Nyquist frequency with respect to the variations of the signal as it reaches the "camera".  Nonetheless, we know that the more points we add to the bandlimited frequency representation, the more accurately we will be able to determine the signal in the space domain.

Once we have finished mapping all points from all n into frequency space, we use an inverse Fourier Transform to obtain our higher resolution image.

On to 2D Tsai and Huang