Re: date-compensated fourier transform

[Rafael Laboissiere <rafael@debian.org>]

* Justin Pryzby <justinpryzby@users.sourceforge.net> [2005-11-17 01:27]:

> Does anyone have experience with date-compensated fourier transforms?
> This is a transform where the input data is spaced unevenly in time.
> I've glanced at some papers, but haven't seen any library code that
> implements this.  I'm interested in finding such code, but it would
> also be nice to understand the algorithm in intuitive terms.

It happens that I was just needing to do this the last week and I came
out with a straightforward way to do it using Octave.  I will try to give
you a tutorial on the method, so launch octave and cut and paste the
text below into your octave prompt.

# A FTT can be seen as a least-square regression.  To convince
# yourself, let us consider the following data:

    x = [10 : -1 : 1]'

# The FFT of x is:

    f = fft (x)
    
# You could obtain the same result using a least square regression:

    A = fft (eye (10)) / 10 ;
    f1 = A \ x
    
# The matrix A contains the cosine/sine basis functions against which the
# regression is done.  To visualize them, try the following (this is for
# a n=100 points FFT, the plot below is only for the first 10 cosine
# components):

    plot (real (fft (eye (100)) (:, 1 : 10)))
    
# Now, imagine that the data has "holes" in it, in other words, we have
# only the points for the following index matrix:

    idx = [1 : 3, 6 : 10];
    plot (idx, x (idx), "*");

# You can still apply the least-square regression by dropping lines in the
# A matrix:

    f2 = A (idx, :) \ x (idx)

# As we expect, f1 and f2 will be different in general.  

# You can see how good both FFT interpolations are by first interpolating
# the A matrix:

    B = fft (eye (100)) (:, [1 : 5, 96 : 100]) / 10;

# and then applying it to both FFT coefficients vectors obtained before:

    x1 = real (B * f1);
    x2 = real (B * f2);
    
# Let us plot it:

    t = linspace (1, 11, 100);
    plot (x, "+r;complete data;", \
          idx, x (idx), "*r;data with holes;", \
	  t, x1, "b;FFT interpolation (complete data);", \
	  t, x2, "c;FFT interpolation (data with holes);");

Hope this is helpful.

-- 
Rafael Laboissiere