Thread Subject: DFT of an irregular time-spacing signal

Hi there!

I have a signal, y which has irregular time-spacing. I would like to find
the power spectrum of this signal. Since the signal has irregular
time-spacing, I guess I can't simply perform the DFT using the fft command
in Matlab?? If so, what method should I used in order to find the power
spectrum of y?

Thanks.
CCC

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message news:<ap102a$6ce$1@scotsman.ed.ac.uk>...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?

help dft

Greg

Hi Greg,

It seems like no such a function in Matlab(dft)??

CCC


"Greg Heath" <heath@ll.mit.edu> wrote in message
news:9e910d7a.0210211421.4fcf9eb@posting.google.com...
> "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
news:<ap102a$6ce$1@scotsman.ed.ac.uk>...
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to
find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft
command
> > in Matlab?? If so, what method should I used in order to find the power
> > spectrum of y?
>
> help dft
>
> Greg

In article <ap23dq$94e$1@scotsman.ed.ac.uk>, Chia.Chong@ee.ed.ac.uk
says...
> Hi Greg,
>
> It seems like no such a function in Matlab(dft)??
>

The Signal Processing Toolbox has the function DFTMTX, a way to
calculate the DFT via a matrix multiplication. This might be way too
memory intensive for you.

--Loren

[Note to Peter Boettcher: Hello Peter... This might qualify as a frequently
asked question and I don't see it on the current list. My answer might be
part of a good answer but it is incomplete. - Ken]

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
news:ap102a$6ce$1@scotsman.ed.ac.uk...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC
>
>
>

Hello Chia,

This thread seems to have headed off into a digression. I'll try to bring it
back on track...

The easiest way to find the spectrum of irregularly sampled data is to
resample it to uniform samples. You can do this with Matlab's interp1
function. The accuracy of your spectrum will depend on the accuracy of the
interpolation. You will want to experiment with several of the interpolation
methods that are available in interp1. I believe that for interpolation with
a limited window (i.e. interpolating a sample value from N nearest
neighbors), the Lagrange interpolation is optimal, but Lagrange is not one
of the choices in interp1 (hint, hint Mathworks).

If interpolation doesn't work there are other schemes available. The
Lomb-Scargle periodogram is often mentioned in relation to this question and
may be more appropriate if your data has very uneven spacing (e.g. very
large or very small spacings). I know that this algorithm is listed in
Numerical Recipes, but I don't have a good on-line reference (with Matlab
code) to point you to. Perhaps someone else following this thread has this
information.

In general, the problem is that the spacing between points determines the
"importance" of a particular point. For example, if several points are very
close together, then small amounts of noise on those measurements will tend
to have a greater effect on inferring the slope of the function (and with
it, the high frequency energy) than the same amounts of noise on
measurements that are further apart.

HTH,

Ken

Ken Davis wrote:
> [Note to Peter Boettcher: Hello Peter... This might qualify as a frequently
> asked question and I don't see it on the current list. My answer might be
> part of a good answer but it is incomplete. - Ken]

I see on the "Books" page that there are a huge number of
Signal Processor books listed. Is this topic covered in
any of them?
http://www.mathworks.com/support/books/index.jsp?category=11

          - Randy

Hi CCC,

If the other suggested alternatives do not help you (such as first
interpolating the data to make it uniformly spaced before using the DFT
function), then you may need to write your own MATLAB function to compute
the DFT at only the irregularly-spaced points that you are interested in.
Basically, you would need to find a definition of DFT in a textbook and
implement it directly as a function in M.

In addition, I will let me colleagues here know about this need for a new
DFT function for irregularly-spaced data.

Thanks,
Steve
The MathWorks

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
news:ap102a$6ce$1@scotsman.ed.ac.uk...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC
>
>
>

Hello,

I believe that there is no simple way to write a DFT for unevenly spaced
data. The DFT/FFT (they are the same function, but the FFT is a fast
algorithm for computing the DFT), requires evenly spaced data in it's
underlying assumptions. If the set of samples is not not evenly spaced,
then, in general, the spectrum will not be a finite, discrete set of
frequencies.

If you take the "naive" approach and just take the inner product of your
sample vector with a collection of complex sinusoids sampled with the same
uneven spacing, you will get a set of numbers, but those numbers will not be
a very good representation of the power spectrum since the unevenly sampled
sinusoids will not form an orthonormal set.

HTH,

Ken

"Steve Conahan" <sconahan@mathworks.com> wrote in message
news:ap3otm$sh2$1@news.mathworks.com...
> Hi CCC,
>
> If the other suggested alternatives do not help you (such as first
> interpolating the data to make it uniformly spaced before using the DFT
> function), then you may need to write your own MATLAB function to compute
> the DFT at only the irregularly-spaced points that you are interested in.
> Basically, you would need to find a definition of DFT in a textbook and
> implement it directly as a function in M.
>
> In addition, I will let me colleagues here know about this need for a new
> DFT function for irregularly-spaced data.
>
> Thanks,
> Steve
> The MathWorks
>
> "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to
find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft
command
> > in Matlab?? If so, what method should I used in order to find the power
> > spectrum of y?
> >
> > Thanks.
> > CCC
> >
> >
> >
>
>

Chia C Chong wrote:
>
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC

Try "nonuniform DFT" in google.

Ken Davis wrote:
> Hello,
>
> I believe that there is no simple way to write a DFT for unevenly spaced
> data. The DFT/FFT (they are the same function, but the FFT is a fast
> algorithm for computing the DFT), requires evenly spaced data in it's
> underlying assumptions. If the set of samples is not not evenly spaced,
> then, in general, the spectrum will not be a finite, discrete set of
> frequencies.
>
> If you take the "naive" approach and just take the inner product of your
> sample vector with a collection of complex sinusoids sampled with the same
> uneven spacing, you will get a set of numbers, but those numbers will not be
> a very good representation of the power spectrum since the unevenly sampled
> sinusoids will not form an orthonormal set.

However, if you properly approximate the INTEGRAL of
exp(j*w*t) * f(t) where f(t) is your function, then
you should get a reasonable estimate of the Fourier
component at w.

Fourier theory says the frequency spacing of the DFT
is governed by the limits of the integral, the time
interval [A,B] over which you sample. w will be
integer multiples of 2*pi/(B-A).

Approximating the integral suggests that you
should calculate something roughly like
exp(j*w*t).*f(t).*diff(t). (diff(t) is
the wrong length so this is not exactly
correct).

          - Randy

From one point of view, the DFT (or FFT) is just a mathematical
transformation. You can apply it on any vector of samples and it will result
in another vector of samples. The question is how to interpret the result of
applying the transformation on a non-uniform sampled signal. There are many
multiresolution spectral analysis systems based on calculating the Fourier
transformation of non-uniform sampled signals. If you would like to go into
such systems you might want to start with the following papers and the
references in them:

1) S.K. Mitra, S. Chakrabarti, and E. Abreu, "Nonuniform Discrete Fourier
Transform and its application in Signal Processing", Proc. EUSIPCO '92, Vol.
2, pp. 909-912.

2) J. Garas, and P. Sommen, "Warped Linear Time Invariant Systems and Their
Application in Audio Signal Processing", ICASSP '99. You can get this paper
in postscript here
http://www.dspalgorithms.com/technology/icassp99.html
or in PDF here
http://www.telecom.tuc.gr/paperdb/icassp99/PDF/SCAN/IC991209.PDF

Regards,

John Garas

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> schreef in bericht
news:ap102a$6ce$1@scotsman.ed.ac.uk...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC
>
>
>

Ken Davis wrote:
>
> Hello,
>
> I believe that there is no simple way to write a DFT for unevenly spaced
> data. The DFT/FFT (they are the same function, but the FFT is a fast
> algorithm for computing the DFT), requires evenly spaced data in it's
> underlying assumptions. If the set of samples is not not evenly spaced,
> then, in general, the spectrum will not be a finite, discrete set of
> frequencies.
>
> If you take the "naive" approach and just take the inner product of your
> sample vector with a collection of complex sinusoids sampled with the same
> uneven spacing, you will get a set of numbers, but those numbers will not be
> a very good representation of the power spectrum since the unevenly sampled
> sinusoids will not form an orthonormal set.
>
> HTH,
>
> Ken
>
Amen! A simple and stringent test of such a "transform" is the ability
to recreate the original set of points. One could devise a figure of
merit for various interpolation methods that uses closeness to exact
reproduction as its measure. Mean square error seems appropriate.

Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Hello,

Approximating the integral in this way is roughly equivalent to
interpolating using a zero-order hold as your interpolation function. You
can probably do better with a better interpolator.

Ken

"Randy Poe" <rpoe@nospam.com> wrote in message
news:ap3sdg0a4o@enews3.newsguy.com...
> Ken Davis wrote:
> > Hello,
> >
> > I believe that there is no simple way to write a DFT for unevenly spaced
> > data. The DFT/FFT (they are the same function, but the FFT is a fast
> > algorithm for computing the DFT), requires evenly spaced data in it's
> > underlying assumptions. If the set of samples is not not evenly spaced,
> > then, in general, the spectrum will not be a finite, discrete set of
> > frequencies.
> >
> > If you take the "naive" approach and just take the inner product of your
> > sample vector with a collection of complex sinusoids sampled with the
same
> > uneven spacing, you will get a set of numbers, but those numbers will
not be
> > a very good representation of the power spectrum since the unevenly
sampled
> > sinusoids will not form an orthonormal set.
>
> However, if you properly approximate the INTEGRAL of
> exp(j*w*t) * f(t) where f(t) is your function, then
> you should get a reasonable estimate of the Fourier
> component at w.
>
> Fourier theory says the frequency spacing of the DFT
> is governed by the limits of the integral, the time
> interval [A,B] over which you sample. w will be
> integer multiples of 2*pi/(B-A).
>
> Approximating the integral suggests that you
> should calculate something roughly like
> exp(j*w*t).*f(t).*diff(t). (diff(t) is
> the wrong length so this is not exactly
> correct).
>
> - Randy
>

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message news:<ap23dq$94e$1@scotsman.ed.ac.uk>...
> Hi Greg,
>
> It seems like no such a function in Matlab(dft)??

Sorry Chia,

The dft.m on our local network was written by a former employee.
The code looks like a straightforward implementation of the equations
given in "help fft".

There are dft algorithms around (in Fortran I think) that are fast.
I have a list somewhere that I got from a recent google search.
I think it led me to some old newsgroup posts on either comp.dsp or
sci.math.num-analysis.

I will relay the info if I can find it(my office mate got tired of my
mess and "reorganized" the office).

I think I used the word "nonequidistant" as a search keyword.

Hope this helps.

Greg

> "Greg Heath" <heath@ll.mit.edu> wrote in message
> news:9e910d7a.0210211421.4fcf9eb@posting.google.com...
> > "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> news:<ap102a$6ce$1@scotsman.ed.ac.uk>...
> > > Hi there!
> > >
> > > I have a signal, y which has irregular time-spacing. I would like to
> find
> > > the power spectrum of this signal. Since the signal has irregular
> > > time-spacing, I guess I can't simply perform the DFT using the fft
> command
> > > in Matlab?? If so, what method should I used in order to find the power
> > > spectrum of y?
> >
> > help dft

heath@ll.mit.edu (Greg Heath) wrote in message news:<9e910d7a.0210211421.4fcf9eb@posting.google.com>...
> "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message news:<ap102a$6ce$1@scotsman.ed.ac.uk>...
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft command
> > in Matlab?? If so, what method should I used in order to find the power
> > spectrum of y?
>
> help dft
>
> Greg

I just replied to this 1st reply but cannot reply to my 2nd reply
directly (;>).

Third reply: As indicated in my previous post, I searched *within
newsgoups* on google using "fourier nonequispaced". I found the
following threads containing references:

sci.math.num-analysis

Non-equidistant discrete data: Fourier transform, Bruce Johnson,
May 28,1999.
FFT with non equidistant sample, Peter Spelluci, Apr 11, 2001
FFT algorithm for "sparse" input, Sean Moore, Mar 7, 1995

Hope this helps,

Greg

Stan Pawlukiewicz <stanp@mitre.org> wrote in message
news:3DB57AE3.3886D935@mitre.org...
> Chia C Chong wrote:
> >
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to
find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft
command
> > in Matlab?? If so, what method should I used in order to find the
power
> > spectrum of y?
> >
> > Thanks.
> > CCC
>
> Try "nonuniform DFT" in google.


Try "lomb scargle periodogram" on google.

"Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
news:ap102a$6ce$1@scotsman.ed.ac.uk...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC

To cut to the chase, as they say, here is the code.
Cheers! -Aj

function X=dft(t,x,f)
% function X=dft(t,x,f)
% Compute DFT (Discrete Fourier Transform) at frequencies given
% in f, given samples x taken at times t:
% X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
% k

t = t(:); % Format 't' into a column vector
x = x(:); % Format 'x' into a column vector
f = f(:); % Format 'f' into a column vector

% It's just this simple:
W = exp( 2*pi*j * f*t');
X = W * x;

Dear all,

Thanks everyone for the references & suggestions.


Regards,
CCC

"AJ "no z" Johnson" <aj.jozhnson@lmco.com> wrote in message
news:ap5rco$4952@cui1.lmms.lmco.com...
> "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to
find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft
command
> > in Matlab?? If so, what method should I used in order to find the power
> > spectrum of y?
> >
> > Thanks.
> > CCC
>
> To cut to the chase, as they say, here is the code.
> Cheers! -Aj
>
> function X=dft(t,x,f)
> % function X=dft(t,x,f)
> % Compute DFT (Discrete Fourier Transform) at frequencies given
> % in f, given samples x taken at times t:
> % X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
> % k
>
> t = t(:); % Format 't' into a column vector
> x = x(:); % Format 'x' into a column vector
> f = f(:); % Format 'f' into a column vector
>
> % It's just this simple:
> W = exp( 2*pi*j * f*t');
> X = W * x;
>
>
>

AJ \"no z\" Johnson wrote:
>
> "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > Hi there!
> >
> > I have a signal, y which has irregular time-spacing. I would like to find
> > the power spectrum of this signal. Since the signal has irregular
> > time-spacing, I guess I can't simply perform the DFT using the fft command
> > in Matlab?? If so, what method should I used in order to find the power
> > spectrum of y?
> >
> > Thanks.
> > CCC
>
> To cut to the chase, as they say, here is the code.
> Cheers! -Aj
>
> function X=dft(t,x,f)
> % function X=dft(t,x,f)
> % Compute DFT (Discrete Fourier Transform) at frequencies given
> % in f, given samples x taken at times t:
> % X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
> % k
>
> t = t(:); % Format 't' into a column vector
> x = x(:); % Format 'x' into a column vector
> f = f(:); % Format 'f' into a column vector
>
> % It's just this simple:
> W = exp( 2*pi*j * f*t');
> X = W * x;

How closely does an evaluation of the spectrum information reproduce the
original data?

Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

"Dale B. Dalrymple" wrote:
> >
> > Try "nonuniform DFT" in google.
>
> Try "lomb scargle periodogram" on google.

I thought you were being a smart ass and had to try it and see.
I'll be damned! :-)


Bob
--

"Things should be described as simply as possible, but no
simpler."

                                             A. Einstein


////////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

 To contribute your unused processor cycles to the fight against
cancer:

     http://www.intel.com/cure

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////////////

Bob Cain wrote:
>
> "Dale B. Dalrymple" wrote:
> > >
> > > Try "nonuniform DFT" in google.
> >
> > Try "lomb scargle periodogram" on google.
>
> I thought you were being a smart ass and had to try it and see.
> I'll be damned! :-)
>
> Bob
> --
>
> "Things should be described as simply as possible, but no
> simpler."
>
> A. Einstein
>
>
Me too. As the chop from Brooklyn said, "Live and loin."

Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Jerry Avins <jya@ieee.org> wrote in message news:<3DB6B41C.51397F4F@ieee.org>...
> AJ \"no z\" Johnson wrote:
> >
> > "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> > news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > Hi there!
> > >
> > > I have a signal, y which has irregular time-spacing. I would like to find
> > > the power spectrum of this signal. Since the signal has irregular
> > > time-spacing, I guess I can't simply perform the DFT using the fft command
> > > in Matlab?? If so, what method should I used in order to find the power
> > > spectrum of y?
> > >
> > > Thanks.
> > > CCC
> >
> > To cut to the chase, as they say, here is the code.
> > Cheers! -Aj
> >
> > function X=dft(t,x,f)
> > % function X=dft(t,x,f)
> > % Compute DFT (Discrete Fourier Transform) at frequencies given
> > % in f, given samples x taken at times t:
> > % X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
> > % k
> >
> > t = t(:); % Format 't' into a column vector
> > x = x(:); % Format 'x' into a column vector
> > f = f(:); % Format 'f' into a column vector
> >
> > % It's just this simple:
> > W = exp( 2*pi*j * f*t');
> > X = W * x;
>
> How closely does an evaluation of the spectrum information reproduce the
> original data?

W = exp(-2*pi*j * f*t'); % Note the minus sign

1. For uniformly spaced f and t

N = length(X)
x = (1/N)* W' * X;
 
2. For arbitrarily spaced f and t

x = W / X;

Using both inverse forms, I compared fft/dft and ifft/idft on

x = t*exp(-t)

and

X = 1./(1+2*pi*j*f).^2

for uniformly spaced t and f.

Plotting indicated no difference.

Hope this helps.

Greg

Jerry Avins <jya@ieee.org> wrote in message
news:<3DB6B41C.51397F4F@ieee.org>...
> AJ \"no z\" Johnson wrote:
> >
> > "Chia C Chong" <Chia.Chong@ee.ed.ac.uk> wrote in message
> > news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > Hi there!
> > >
> > > I have a signal, y which has irregular time-spacing. I would like to find
> > > the power spectrum of this signal. Since the signal has irregular
> > > time-spacing, I guess I can't simply perform the DFT using the fft command
> > > in Matlab?? If so, what method should I used in order to find the power
> > > spectrum of y?
> > >
> > > Thanks.
> > > CCC
> >
> > To cut to the chase, as they say, here is the code.
> > Cheers! -Aj
> >
> > function X=dft(t,x,f)
> > % function X=dft(t,x,f)
> > % Compute DFT (Discrete Fourier Transform) at frequencies given
> > % in f, given samples x taken at times t:
> > % X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
> > % k
> >
> > t = t(:); % Format 't' into a column vector
> > x = x(:); % Format 'x' into a column vector
> > f = f(:); % Format 'f' into a column vector
> >
> > % It's just this simple:
> > W = exp( 2*pi*j * f*t');

To conform to the MATLAB convention

 W = exp(-2*pi*j * f*t'); % i.e., notice the minus sign

> > X = W * x;
>
> How closely does an evaluation of the spectrum information reproduce
> the original data?

Obviously, when t and f are uniformly sampled,

N = length(X)
x = (1/N)* W' * X;

otherwise

x = W \ X;

W should be well conditioned as long as min(dt) and min(df) are
not too close to 0. Maybe the limit is given by something like

 min(dt)* min(df) >= 1/N?

Other stuff comes to mind:

1. Is the Nyquist frequency FNy = 1/(2*min(dt))?
2. If so, how does that constrain our choice of f?
    a. max(f) > 2*FNy = 1/min(dt)?
    b. min(df) < 1/max(t+dt)?
3. The series representation of x isn't periodic.
4. How is the 0 centered spectrum obtained?(Can't use fftshift)
5. For a given vector dt, is there an analytic formula for an
   optimal choice of df so that the inverse of W has a simple form?
6. Is

    X = dt*fft(x); x = (1/dt)*real(ifft(X));

generalized to

    X = dft(x.*dt); x = N*real(idft(X.*df)); ?

Hope that helps.

Greg

"Greg Heath" <heath@ll.mit.edu> wrote in message
news:9e910d7a.0210232227.78ca98ff@posting.google.com...
> > > W = exp( 2*pi*j * f*t');
>
> To conform to the MATLAB convention
>
> W = exp(-2*pi*j * f*t'); % i.e., notice the minus sign
>
[snip]
> Greg

Good catch... it's been a while since I used this. I don't think the minus
sign is merely a MATLAB convention. I suspect I may have taken out the minus
sign to perform a crude IDFT, neglecting to put it back in.
-Aj

Greg Heath wrote:
>
> Jerry Avins <jya@ieee.org> wrote in message
  ...
> >
> > How closely does an evaluation of the spectrum information reproduce
> > the original data?
>
> Obviously, when t and f are uniformly sampled,
>
> N = length(X)
> x = (1/N)* W' * X;
>
> otherwise
>
> x = W \ X;
>
> W should be well conditioned as long as min(dt) and min(df) are
> not too close to 0. Maybe the limit is given by something like
>
> min(dt)* min(df) >= 1/N?
>
> Other stuff comes to mind:
>
> 1. Is the Nyquist frequency FNy = 1/(2*min(dt))?
> 2. If so, how does that constrain our choice of f?
> a. max(f) > 2*FNy = 1/min(dt)?
> b. min(df) < 1/max(t+dt)?
> 3. The series representation of x isn't periodic.
> 4. How is the 0 centered spectrum obtained?(Can't use fftshift)
> 5. For a given vector dt, is there an analytic formula for an
> optimal choice of df so that the inverse of W has a simple form?
> 6. Is
>
> X = dt*fft(x); x = (1/dt)*real(ifft(X));
>
> generalized to
>
> X = dft(x.*dt); x = N*real(idft(X.*df)); ?
>
> Hope that helps.
>
> Greg

The sine and cosine coefficients from the analysis (any form of the
result can be recast that way) define a continuous function that is
claimed to represent the signal that was sampled. When the original
signal is bandlimited and uniformly sampled at a high enough rate, that
function indeed represents the actual samples when evaluated at the
sampling instants. If the analysis for non-uniform sampling is valid,
that condition still must hold. It is not too difficult to imagine
sample sets in which it doesn't hold. The interesting problem is do find
the properties of sets for which it does. I suspect without proof that
requiring the longest sampling interval to be no longer than the period
needed for uniform sampling is sufficient, but unnecessarily
restrictive.

Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Jerry Avins wrote:
>
  ...
>
> The sine and cosine coefficients from the analysis (any form of the
> result can be recast that way) define a continuous function that is
> claimed to represent the signal that was sampled. When the original
> signal is bandlimited and uniformly sampled at a high enough rate, that
> function indeed represents the actual samples when evaluated at the
> sampling instants. If the analysis for non-uniform sampling is valid,
> that condition still must hold. It is not too difficult to imagine
> sample sets in which it doesn't hold. The interesting problem is do find
> the properties of sets for which it does. I suspect without proof that
> requiring the longest sampling interval to be no longer than the period
> needed for uniform sampling is sufficient, but unnecessarily
> restrictive.
>
> Jerry
> --
... and got an email response. My response to that is reproduced without
revealing the writer:

Your argument is valid (I think) with complete knowledge. Given noise,
measurement uncertainty, and quantization, it becomes increasingly worse
as some measurements bunch together. If no sampling interval exceeds an
acceptable interval for uniform sampling, there must be a good solution.
samples taken much closer in time that the period of the highest
frequency won't be completely independent, so more samples are needed
than the minimum with uniform sampling. Extreme case: Signal band
limited to 5 KHz; 20,000 samples taken over 1 second. With uniform
sampling, this is 2X oversampled. If the samples are uniformly
distributed over the first 2 milliseconds, leaving the last 998
milliseconds as a sort of terra incognita, reconstruction is difficult
at best. If we replace "milliseconds" with "microseconds", the
mathematical precision needed for reconstruction is unobtainable. More
comments below.

Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
********** wrote:
>
> >
> > The sine and cosine coefficients from the analysis (any form of the
> > result can be recast that way) define a continuous function that is
> > claimed to represent the signal that was sampled. When the original
> > signal is bandlimited and uniformly sampled at a high enough rate, that
> > function indeed represents the actual samples when evaluated at the
> > sampling instants. If the analysis for non-uniform sampling is valid,
> > that condition still must hold. It is not too difficult to imagine
> > sample sets in which it doesn't hold. The interesting problem is do find
> > the properties of sets for which it does. I suspect without proof that
> > requiring the longest sampling interval to be no longer than the period
> > needed for uniform sampling is sufficient, but unnecessarily
> > restrictive.
> >
> > Jerry
> > --
> > Engineering is the art of making what you want from things you can get.
> > ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
>
> (good quote)

I'm glad you like it. It's original.
 
> Here are some thoughts, without actually doing any analysis.
> Think of the DFT as solving a set of simulaneous equations. The time-based
> data points are the constants, and we are solving for the coefficients that
> represent the magnitudes of certain sinusoids. If there are N data points,
> there are N*2 constants (including the imaginary parts), and we are solving
> for N*2 unknowns (N coefficients for each of sine and cosine). Assuming the
> set of sinusoids is a linearly independent, the solution is unique.
>
> In addition, because of the assumption that the input data was band-limited,
> the DFT coefficients are guaranteed to reproduce the original signal,
> exactly.
>
> But now say we have N data point that are spaced irregularly. The set of
> linear equations still has a unique solution (?). The question is, to what
> extent can the original waveform be reproduced?
>
> Certainly, the inverse DFT will produce a waveform that passes through each
> of the input points, and the inverse will also be bandlimited, since it is
> composed of a limear combination of band-limited sines and cosines. So if
> the input was bandlimited, the original signal is exactly reconstructed.
>
> This argument has a flaw; I think there are cases where the simultaneous
> equations do not have full rank. Say, for example, the data is sampled at
> exactly twice the normal interval. Well then, the DFT implicity "assumes"
> the missing points are zero; but really, there is a sinusoid at some
> harmonic frequency whose magnitude is indeterminate!

But DFT probably doesn't work at all. N equations, N unknowns: actual
solution method undefined.
 
> I have a hunch that the signal reconstruction is determined more by the
> shortest interval, rather than the longest, or perhaps it is sufficient that
> the intervals have some sort of relatively-prime relationship.

You mean that if two samples are taken sufficiently close, a useful
solution is guaranteed? Surely not!
 
Jerry
--
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Jerry Avins <jya@ieee.org> wrote in message
news:<3DB84626.4A02055C@ieee.org>...
> Greg Heath wrote:
> >
> > Jerry Avins <jya@ieee.org> wrote in message
> ...
> > >
> > > How closely does an evaluation of the spectrum information
> > > reproduce the original data?
> >
> > Obviously, when t and f are uniformly sampled,
> >
> > N = length(X)
> > x = (1/N)* W' * X;
> >
> > otherwise
> >
> > x = W \ X;
> >
> > W should be well conditioned as long as min(dt) and min(df) are
> > not too close to 0. Maybe the limit is given by something like
> >
> > min(dt)* min(df) >= 1/N?
> >
> > Other stuff comes to mind:
> >
> > 1. Is the Nyquist frequency FNy = 1/(2*min(dt))?
> > 2. If so, how does that constrain our choice of f?
> > a. max(f) > 2*FNy = 1/min(dt)?
> > b. min(df) < 1/max(t+dt)?
> > 3. The series representation of x isn't periodic.
> > 4. How is the 0 centered spectrum obtained?(Can't use fftshift)
> > 5. For a given vector dt, is there an analytic formula for an
> > optimal choice of df so that the inverse of W has a simple form?
> > 6. Is
> >
> > X = dt*fft(x); x = (1/dt)*real(ifft(X));
> >
> > generalized to
> >
> > X = dft(x.*dt); x = N*real(idft(X.*df)); ?
>
> The sine and cosine coefficients from the analysis (any form of the
> result can be recast that way) define a continuous function that is
> claimed to represent the signal that was sampled. When the original
> signal is bandlimited and uniformly sampled at a high enough rate, that
> function indeed represents the actual samples when evaluated at the
> sampling instants. If the analysis for non-uniform sampling is valid,
> that condition still must hold. It is not too difficult to imagine
> sample sets in which it doesn't hold. The interesting problem is do find
> the properties of sets for which it does. I suspect without proof that
> requiring the longest sampling interval to be no longer than the period
> needed for uniform sampling is sufficient, but unnecessarily
> restrictive.

If you mean for bandwidth B, max(dt) < 1/B should be sufficient,
then, I agree. Furthermore, since

T = sum(k=1,N){dt(k)} <= N*max(dt),

a frequency resolution lower bound for *uniform* freqency
sampling is

df = B/N < [1/(N*max(dt))] <= 1/T,

which corresponds to 2b (above).


   However, isn't

min(dt) < 1/B

(contradicting 2a above) always necessary? Then if

min(dt) < 1/B < max(dt),

we get

[1/(N*max(dt))]< df = B/N < [1/(N*min(dt))].

Therefore, the choice

df = 1/T

is still allowed.

Hope this helps,

Greg

heath@ll.mit.edu (Greg Heath) wrote in message

news:<9e910d7a.0210251842.4616d1a1@posting.google.com>...
> Jerry Avins <jya@ieee.org> wrote in message
> news:<3DB84626.4A02055C@ieee.org>...
> > Greg Heath wrote:
> > >
> > > Jerry Avins <jya@ieee.org> wrote in message
> > > > ...
> > > > How closely does an evaluation of the spectrum information
> > > > reproduce the original data?
> > >
> > > Obviously, when t and f are uniformly sampled,
> > >
> > > N = length(X)
> > > x = (1/N)* W' * X;
> > >
> > > otherwise
> > >
> > > x = W \ X;
> > >
> > > W should be well conditioned as long as min(dt) and min(df) are
> > > not too close to 0. Maybe the limit is given by something like
> > >
> > > min(dt)* min(df) >= 1/N?
> > >
> > > Other stuff comes to mind:
> > >
> > > 1. Is the Nyquist frequency FNy = 1/(2*min(dt))?

No. In fact, it can be *orders of magnitude* larger (Scargle,1982)!
Let

t = (t1,...,tN),
  = t1 + k*dt0,

where dt0 is the largest common multiple of ti-t1 and k = (k1,...,kN) is the
resulting vector of integers. Then

FNy = 1/(2*dt0) >= 1/(2*min(dt)).

Typically, dt0 << min(dt) (i.e., the ki are very large). Then

FNy = 1/(2*dt0) >> 1/(2*min(dt))!

For example, if min(dt) = 0.01 and the ti are stored to 8-digit accuracy,
dt0 can be as small as 1e-8. Then

FNy = 1/(2*1e-8) = 50 MHz >> 1/(2*0.01) = 50 Hz.

Two important points: With nonuniform sampling,
1. the Nyquist frequency is a function of the sampling pattern and can be
   unrelated to the bandwidth of the underlying continuous time function.
2. The Nyquist frequency can be increased by orders of magnitude by using
   nonuniform sampling.

For example, if the last digit is truncated in the above example,
dt0 = 1e-7 and FNy = 5 MHz.

> > > 2. If so, how does that constrain our choice of f?
> > > a. max(f) > 2*FNy = 1/min(dt)?

To avoid aliasing

max(f) < 2*FNy = 1/(dt0).

Therefore, it is typical, with nonuniform sampling, to have
  
max(f) >> 1/min(dt).

> > > b. min(df) < 1/max(t+dt)?

Choose df to be uniform and define dtN to be a reasonable
(i.e., min(dt) <= dtN <= max(dt)), integer multiple of dt0. Then there is
no loss of frequency resolution if

df <= 1/T

T = tN + dtN = M * dt0

> > > 3. The series representation of x isn't periodic.

Incorrect. If df = 1/(M*dt0), x will have the period T.

> > > 4. How is the 0 centered spectrum obtained?(Can't use fftshift)

By determining dt0 and fNy.

> > > 5. For a given vector dt, is there an analytic formula for an
> > > optimal choice of df so that the inverse of W has a simple form?

Not exactly. A common choice is to zero pad x and use a power of 2 ifft.
An interpolation of the original sampled points will be recovered.

Greg

Hope this helps.

Gregory E. Heath heath@ll.mit.edu The views expressed here are
M.I.T. Lincoln Lab (781) 981-2815 not necessarily shared by
Lexington, MA (781) 981-0908(FAX) M.I.T./LL or its sponsors
02420-9185, USA


> > > 6. Is
> > >
> > > X = dt*fft(x); x = (1/dt)*real(ifft(X));
> > >
> > > generalized to
> > >
> > > X = dft(x.*dt); x = N*real(idft(X.*df)); ?
> >
> > The sine and cosine coefficients from the analysis (any form of the
> > result can be recast that way) define a continuous function that is
> > claimed to represent the signal that was sampled. When the original
> > signal is bandlimited and uniformly sampled at a high enough rate, that
> > function indeed represents the actual samples when evaluated at the
> > sampling instants. If the analysis for non-uniform sampling is valid,
> > that condition still must hold. It is not too difficult to imagine
> > sample sets in which it doesn't hold. The interesting problem is do find
> > the properties of sets for which it does. I suspect without proof that
> > requiring the longest sampling interval to be no longer than the period
> > needed for uniform sampling is sufficient, but unnecessarily
> > restrictive.
>
> If you mean for bandwidth B, max(dt) < 1/B should be sufficient,
> then, I agree. Furthermore, since
>
> T = sum(k=1,N){dt(k)} <= N*max(dt),
>
> a frequency resolution lower bound for *uniform* freqency
> sampling is
>
> df = B/N < [1/(N*max(dt))] <= 1/T,
>
> which corresponds to 2b (above).
>
>
> However, isn't
>
> min(dt) < 1/B
>
> (contradicting 2a above) always necessary? Then if
>
> min(dt) < 1/B < max(dt),
>
> we get
>
> [1/(N*max(dt))]< df = B/N < [1/(N*min(dt))].
>
> Therefore, the choice
>
> df = 1/T
>
> is still allowed.
>
> Hope this helps,
>
> Greg

What do you mean by "irregular" time-spacing? If you can characterize the
statistics of the sampling function then you can also calculate the PSD of
the sampled signal directly, using stochastic methods.

The relations between frequency and time are (DFT, FFT, etc.) are a
convenient tool for many applications, but not the only tool in the box.
For example, you could use time-series estimation techniques with selective
mapping of samples to sampling bins.

This is generally discussed in books on multi-rate sampling systems. I
would need more details regarding the application to be of more help.

- Michael Pender

Executive Manager
Nanochron, LLC


Chia C Chong <Chia.Chong@ee.ed.ac.uk> wrote in message
news:ap102a$6ce$1@scotsman.ed.ac.uk...
> Hi there!
>
> I have a signal, y which has irregular time-spacing. I would like to find
> the power spectrum of this signal. Since the signal has irregular
> time-spacing, I guess I can't simply perform the DFT using the fft command
> in Matlab?? If so, what method should I used in order to find the power
> spectrum of y?
>
> Thanks.
> CCC
>
>
>

In article <1R3w9.38629$iV1.8275@nwrddc02.gnilink.net>,
 "Michael Pender" <michael.pender@nanochron.com> writes:
>What do you mean by "irregular" time-spacing? If you can characterize the

snip

>Chia C Chong <Chia.Chong@ee.ed.ac.uk> wrote in message
>news:ap102a$6ce$1@scotsman.ed.ac.uk...
>> Hi there!
>>
>> I have a signal, y which has irregular time-spacing. I would like to find
>> the power spectrum of this signal. Since the signal has irregular
>> time-spacing, I guess I can't simply perform the DFT using the fft command
>> in Matlab?? If so, what method should I used in order to find the power
>> spectrum of y?
>>

snip

from my annotations:
\begin{citation}
You might take a look at "Spectral Analysis of Unevenly Sampled Data"
in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
chapter will instruct you in using a periodogram-approach (called the
Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
welcome to download it and give it a whirl.
Regards,
Brett Shoelson
\end{citation}

\begin{citation}
Luis <lfarina@uvigo.es> wrote:
> Can I do a Fourier Transform of data thats not equaly spaced in time?
> If so, is it the same as the normal Fourier Transform or do I have to
> account for anything more?

There are some articles published by A. Dutt an V. Rokhlin which deal
with this problem. For a start you can try

\bibitem{dutt1}
A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.

\bibitem{dutt2}
A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
85--100.

These articles contain algorithms for different types of FTs with non-
equispaced data. The implementation of these algorithms isn't trivial,
but they do work.
\end{citation}

\begin{citation}
we are pleased to announce the availability of Version 1.0 of a
C library for computing the Nonequispaced Discrete Fourier Transform
(NDFT) in one, two or three dimensions.
Other common names for NFFT are non-uniform fast Fourier transform
(NUFFT), generalized fast Fourier transform (GFFT),
unequally-spaced fast Fourier transform (USFFT),
non-equispaced fast Fourier transform (NFFT), or gridding.

Our library is free software based on FFTW.

Visit our NFFT web-page

         http://www.math.uni-luebeck.de/potts/nfft

for software, documentation, and related links.


For further information please contact
Daniel Potts (potts@math.uni-luebeck.de).
\end{citation}

hth
peter

On Mon, 21 Oct 2002 14:39:45 +0100, "Chia C Chong"
<Chia.Chong@ee.ed.ac.uk> wrote:

>Hi there!
>
>I have a signal, y which has irregular time-spacing. I would like to find
>the power spectrum of this signal. Since the signal has irregular
>time-spacing, I guess I can't simply perform the DFT using the fft command
>in Matlab?? If so, what method should I used in order to find the power
>spectrum of y?

Maybe this could be a solution. Try ;-)

%F_LINFOURIER Harmonic analysis using a least square method
%
% [R,phi,A,B,yr]=f_linfourier(t,y,f);
%
% t: time
% y: signal
% f: array of frequencies
% A: amplitude of cosinusoidal components
% B: amplitude of sinusoidal components
% R: (A.^2+B.^2).^0.5 but if f(k)=0, then R(k)=A(k)=mean(y);
% phi: atan2(-B,A) phase relative to the initial time
% yr: reconstructed signal
%
% The approximate signal at time t could be computed in this way:
% yr=NaN*zeros(length(t),2*length(f));
%
f_mat=repmat(f(:).',length(t),2);t_rel=t-t(1);t_mat=repmat(t_rel(:),1,2*length(f));
% M=NaN*zeros(length(t),2*length(f));
%
M(:,1:length(f))=cos(2*pi*f_mat(:,1:length(f)).*t_mat(:,1:length(f)));
%
M(:,length(f)+1:end)=sin(2*pi*f_mat(:,length(f)+1:end).*t_mat(:,length(f)+1:end));
% yr=sum(M*[A(:),B(:)],2);
%
% the previous lines are the code for the following formula:
%
% N
% yr= sum A(k)*cos(2*pi*f(k)*(t-t(1)))+B(k)*sin(2*pi*f(k)*(t-t(1)))
% k=1
%
% or, it is equivalent,
%
% N
% yr= sum R(k)*cos(2*pi*f(k)*(t-t(1))+phi(k))
% k=1
%
% Where N=length(f)
%
% See also:
% FFT, F_FFTFREQ, F_FFTMOB
%
% Gabriele Bulian (ruga@libero.it) 24/04/2002

function [R,phi,A,B,yr]=f_linfourier(t,y,f);

t=(t-t(1));

f=f(:).';
y=y(:);
t=t(:);
N=length(t);
Nf=length(f);
f_mat=repmat(f,N,2);
t_mat=repmat(t,1,2*Nf);
M=NaN*zeros(N,2*Nf);
M(:,1:Nf)=cos(2*pi*f_mat(:,1:Nf).*t_mat(:,1:Nf));
M(:,Nf+1:end)=sin(2*pi*f_mat(:,Nf+1:end).*t_mat(:,Nf+1:end));
S=M\y;
A=S(1:Nf);
B=S(Nf+1:end);

R=(A.^2+B.^2).^0.5;
phi=atan2(-B,A);

R(f==0)=R(f==0).*cos(phi(f==0)); %take into account the sign of A(k)
directly in R
phi(f==0)=0; %and set the phase to zero (otherwise the phase could be
0 or pi or -pi)

if nargout==5,
   yr=sum(M*S,2);
end;




Ciao e 73-51 de Tartaruga .

.oO-=> TARTARUGA (* Gabriele *) <=-Oo.

E-Mail: ruga@libero.it

Packet CB: http://www.geocities.com/ita490/
Astronomy: http://www.geocities.com/rugabri/

"Chi dorme non piglia pesci, ma chi non dorme, alla fin fine...muore..."
(C) Tartaruga 1999 ;-)

"Michael Pender" <michael.pender@nanochron.com> wrote in message news:<1R3w9.38629$iV1.8275@nwrddc02.gnilink.net>...
> What do you mean by "irregular" time-spacing? If you can characterize the
> statistics of the sampling function then you can also calculate the PSD of
> the sampled signal directly, using stochastic methods.

In my case the irregularity cannot be quantified: Radar pulse spacing
conditional on the number, speed and range of targets within the field
of view.

Hope this helps.

Greg

ruga@libero.it (Gabriele) wrote in message news:<3dc4323c.33575065@powernews.libero.it>...
> On Mon, 21 Oct 2002 14:39:45 +0100, "Chia C Chong"
> <Chia.Chong@ee.ed.ac.uk> wrote:
>
> >Hi there!
> >
> >I have a signal, y which has irregular time-spacing. I would like to find
> >the power spectrum of this signal. Since the signal has irregular
> >time-spacing, I guess I can't simply perform the DFT using the fft command
> >in Matlab?? If so, what method should I used in order to find the power
> >spectrum of y?
>
> Maybe this could be a solution. Try ;-)
>
> %F_LINFOURIER Harmonic analysis using a least square method
> %
> % [R,phi,A,B,yr]=f_linfourier(t,y,f);
> %
> % t: time
> % y: signal
> % f: array of frequencies
> % A: amplitude of cosinusoidal components
> % B: amplitude of sinusoidal components
> % R: (A.^2+B.^2).^0.5 but if f(k)=0, then R(k)=A(k)=mean(y);
> % phi: atan2(-B,A) phase relative to the initial time
> % yr: reconstructed signal
> %
> % The approximate signal at time t could be computed in this way:
> % yr=NaN*zeros(length(t),2*length(f));
> %
> f_mat=repmat(f(:).',length(t),2);t_rel=t-t(1);t_mat=repmat(t_rel(:),1,2*length(f));
> % M=NaN*zeros(length(t),2*length(f));
> %
> M(:,1:length(f))=cos(2*pi*f_mat(:,1:length(f)).*t_mat(:,1:length(f)));
> %
> M(:,length(f)+1:end)=sin(2*pi*f_mat(:,length(f)+1:end).*t_mat(:,length(f)+1:end));
> % yr=sum(M*[A(:),B(:)],2);
> %
> % the previous lines are the code for the following formula:
> %
> % N
> % yr= sum A(k)*cos(2*pi*f(k)*(t-t(1)))+B(k)*sin(2*pi*f(k)*(t-t(1)))
> % k=1
> %
> % or, it is equivalent,
> %
> % N
> % yr= sum R(k)*cos(2*pi*f(k)*(t-t(1))+phi(k))
> % k=1
> %
> % Where N=length(f)
> %
> % See also:
> % FFT, F_FFTFREQ, F_FFTMOB
> %
> % Gabriele Bulian (ruga@libero.it) 24/04/2002
>
> function [R,phi,A,B,yr]=f_linfourier(t,y,f);
>
> t=(t-t(1));
>
> f=f(:).';
> y=y(:);
> t=t(:);
> N=length(t);
> Nf=length(f);
> f_mat=repmat(f,N,2);
> t_mat=repmat(t,1,2*Nf);
> M=NaN*zeros(N,2*Nf);
> M(:,1:Nf)=cos(2*pi*f_mat(:,1:Nf).*t_mat(:,1:Nf));
> M(:,Nf+1:end)=sin(2*pi*f_mat(:,Nf+1:end).*t_mat(:,Nf+1:end));
> S=M\y;
> A=S(1:Nf);
> B=S(Nf+1:end);
>
> R=(A.^2+B.^2).^0.5;
> phi=atan2(-B,A);
>
> R(f==0)=R(f==0).*cos(phi(f==0)); %take into account the sign of A(k)
> directly in R
> phi(f==0)=0; %and set the phase to zero (otherwise the phase could be
> 0 or pi or -pi)
>
> if nargout==5,
> yr=sum(M*S,2);
> end;
>
>
>
>
> Ciao e 73-51 de Tartaruga .
>
> .oO-=> TARTARUGA (* Gabriele *) <=-Oo.
>
> E-Mail: ruga@libero.it
>
> Packet CB: http://www.geocities.com/ita490/
> Astronomy: http://www.geocities.com/rugabri/
>
> "Chi dorme non piglia pesci, ma chi non dorme, alla fin fine...muore..."
> (C) Tartaruga 1999 ;-)

Did you see the 23 oct code of AJ Johnson before posting this?

Greg

Here is my 2 cents, lest anyone think there is a solution to the
problem of "spectral analysis of gappy data".


The "ideal" solution to calculating a spectrum of an irregularly
spaced time series would be a least squares fit to the "basis functions"
of sinusoids. In the regularly sampled time series case, the FFT is in fact the
least squares fit solution (with an error = 0, since it is an orthogonal transformation).
If there are gaps, or if the time series is irregularly
sampled, the fourier sampled sinusoids are not longer orthogonal. In fact,
the least squares fit is, in general, ill conditioned. In the dim past, I had spent
a long time trying to find a combination of sinusoids that would reproduce an
orthogonal subspace for "gappy timeseries" but could not find any.


Lomb Scargle is a least squares fit to SINGLE (complex) sinusoid. It is then repeated
to calculated the parameters of other sinusoids. It does not simultaneously
solve the least squares fit to a spectrum. Don't think that lomb periodogram
allows one to calculate the spectrum of a gappy time series. (I had at one time
a piece of IDL code that reproduced the lomb periodogram of a gappy time series
by using an fft (of the same gappy data) by employing the appropriate sampling and
normalization). IMHO the real utility of a lomb aproach is that it can give you
the confidence interval of the amplitude of a harmonic component of an irregularly
sampled time series. IMHO it should not be used to reproduce a power spectrum of an
irregularly sampled time series.

Directly applying the DFT is also inadequate. There is no difference to the case of
putting zeroes in the gaps and FFTing. (In fact, when applying this technique, the user
implicitly puts zeros into the gaps, and perhaps applies different scaling). If you think of
how you calculate the coefficient of a spectral component, it is jst the inner product of
the time series with the sinusoid. If you have gaps, those points are left out of the
inner product. This is identical to placing zeroes in the summation. The only difference
may be the normalization of the inner product, and this does not solve the problem of
gaps, it is merely a different normalzation.

One person shows that this approach, when the time series is calculated by "inverting"
the transform, does indeed produce the time series. That is because they really just
did an FFT of zero-padded gapped data.

Another approach is in calculating the autocorrelation of the time series (which has
no gaps) and FFT ing that to get the power spectrum. This does not work either.



To get a feel of where I am coming from, write out the equation for the spectrum
Ax = b
where b is the time series, A is the transformation matrix, and x is the spectrum. Solve for x.
If you have gappy/irregular data, the columns of A are not orthogonal, and you cannot
immediately solve for x as x = A^T b (like you do in the usual DFT case, and of course, the
FFT is just a fast way of computing A^T b).
Any attempts to force a DFT calculation of gappy data is really just approximating
A^(-1) with A^T which will not, in general, be true. Lomb scargle, which only does one sinusoid
at a time, does this. It assumes A^(-1) = A^T) as well. And as I mentioned before, A^(-1) is
a poorly conditioned problem for typical gappy cases (or at least I have found that in
my particular cases).
What seems to be a good idea is to reduce the rank of A and solve for an underdetermined
case, (i.e. reduce the columns of A to less than the rows of A). I haven't found this to
be a good general solution, but of course that depends on the particulars of the A matrix.


So, in conclusion, the least squares fit (simultaneously for all sinusoids) sounds like
the best approach to me, and all the other approaches (direct DFT calc, and lomb scargle etc)
are inadequate (you might as well just zero padd the gaps and fft).
There is no "good" solution to spectral estimation of irregularly sampled time series.

Cheers,
bob stockwell

On 5 Nov 2002 00:32:26 -0800, heath@ll.mit.edu (Greg Heath) wrote:

>Did you see the 23 oct code of AJ Johnson before posting this?

No, sorry, I was off this group from june to about 1st November.
Should I have seen it ?

>Greg


Ciao e 73-51 de Tartaruga .

.oO-=> TARTARUGA (* Gabriele *) <=-Oo.

E-Mail: ruga@libero.it

Packet CB: http://www.geocities.com/ita490/
Astronomy: http://www.geocities.com/rugabri/

"Chi dorme non piglia pesci, ma chi non dorme, alla fin fine...muore..."
(C) Tartaruga 1999 ;-)

ruga@libero.it (Gabriele) wrote in message
news:<3dcd9ca1.18907218@powernews.libero.it>...
> On 5 Nov 2002 00:32:26 -0800, heath@ll.mit.edu (Greg Heath) wrote:
>
> >Did you see the 23 oct code of AJ Johnson before posting this?
>
> No, sorry, I was off this group from june to about 1st November.
> Should I have seen it ?

Yes(I inserted a minus sign in the exponential):

function X=dft(t,x,f)
% function X=dft(t,x,f)
% Compute DFT (Discrete Fourier Transform) at frequencies given
% in f, given samples x taken at times t:
% X(f) = sum { x(k) * e**(-2*pi*j*t(k)*f) }
% k

t = t(:); % Format 't' into a column vector
x = x(:); % Format 'x' into a column vector
f = f(:); % Format 'f' into a column vector

% It's just this simple:
W = exp(-2*pi*j * f*t');
X = W * x;

Greg

I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.

spellucci@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0i$1@fb04373.mathematik.tu-darmstadt.de>...
>
> In article <1R3w9.38629$iV1.8275@nwrddc02.gnilink.net>,
> "Michael Pender" <michael.pender@nanochron.com> writes:
> >What do you mean by "irregular" time-spacing? If you can characterize the
>
> snip
>
> >Chia C Chong <Chia.Chong@ee.ed.ac.uk> wrote in message
> >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> >> Hi there!
> >>
> >> I have a signal, y which has irregular time-spacing. I would like to find
> >> the power spectrum of this signal. Since the signal has irregular
> >> time-spacing, I guess I can't simply perform the DFT using the fft command
> >> in Matlab?? If so, what method should I used in order to find the power
> >> spectrum of y?
> >>
>
> snip
>
> from my annotations:
> \begin{citation}
> You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> chapter will instruct you in using a periodogram-approach (called the
> Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> welcome to download it and give it a whirl.
> Regards,
> Brett Shoelson
> \end{citation}
>
> \begin{citation}
> Luis <lfarina@uvigo.es> wrote:
> > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > If so, is it the same as the normal Fourier Transform or do I have to
> > account for anything more?
>
> There are some articles published by A. Dutt an V. Rokhlin which deal
> with this problem. For a start you can try
>
> \bibitem{dutt1}
> A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
>
> \bibitem{dutt2}
> A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> 85--100.
>
> These articles contain algorithms for different types of FTs with non-
> equispaced data. The implementation of these algorithms isn't trivial,
> but they do work.
> \end{citation}
>
> \begin{citation}
> we are pleased to announce the availability of Version 1.0 of a
> C library for computing the Nonequispaced Discrete Fourier Transform
> (NDFT) in one, two or three dimensions.
> Other common names for NFFT are non-uniform fast Fourier transform
> (NUFFT), generalized fast Fourier transform (GFFT),
> unequally-spaced fast Fourier transform (USFFT),
> non-equispaced fast Fourier transform (NFFT), or gridding.
>
> Our library is free software based on FFTW.
>
> Visit our NFFT web-page
>
> http://www.math.uni-luebeck.de/potts/nfft
>
> for software, documentation, and related links.
>
>
> For further information please contact
> Daniel Potts (potts@math.uni-luebeck.de).
> \end{citation}
>
> hth
> peter

On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
>
> spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
>
> > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > >What do you mean by "irregular" time-spacing? If you can characterize the
>
> > snip
>
> > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > >> Hi there!
>
> > >> I have a signal, y which has irregular time-spacing. I would like to find
> > >> the power spectrum of this signal. Since the signal has irregular
> > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > >> in Matlab?? If so, what method should I used in order to find the power
> > >> spectrum of y?
>
> > snip
>
> > from my annotations:
> > \begin{citation}
> > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > chapter will instruct you in using a periodogram-approach (called the
> > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > welcome to download it and give it a whirl.
> > Regards,
> > Brett Shoelson
> > \end{citation}
>
> > \begin{citation}
> > Luis <lfar...@uvigo.es> wrote:
> > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > If so, is it the same as the normal Fourier Transform or do I have to
> > > account for anything more?
>
> > There are some articles published by A. Dutt an V. Rokhlin which deal
> > with this problem. For a start you can try
>
> > \bibitem{dutt1}
> > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
>
> > \bibitem{dutt2}
> > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > 85--100.
>
> > These articles contain algorithms for different types of FTs with non-
> > equispaced data. The implementation of these algorithms isn't trivial,
> > but they do work.
> > \end{citation}
>
> > \begin{citation}
> > we are pleased to announce the availability of Version 1.0 of a
> > C library for computing the Nonequispaced Discrete Fourier Transform
> > (NDFT) in one, two or three dimensions.
> > Other common names for NFFT are non-uniform fast Fourier transform
> > (NUFFT), generalized fast Fourier transform (GFFT),
> > unequally-spaced fast Fourier transform (USFFT),
> > non-equispaced fast Fourier transform (NFFT), or gridding.
>
> > Our library is free software based on FFTW.
>
> > Visit our NFFT web-page
>
> > http://www.math.uni-luebeck.de/potts/nfft
>
> > for software, documentation, and related links.
>
> > For further information please contact
> > Daniel Potts (po...@math.uni-luebeck.de).
> > \end{citation}
>
> > hth
> > peter

You might find the more recent discussions on this here in comp.soft-
sys.matlab in April-May of 2008 interesting. Useful is a different
matter.

Dale B. Dalrymple

dbd <dbd@ieee.org> wrote in message <9d1ee7a3-71c6-405f-87a6-0a7fddb46b68@j9g2000prh.googlegroups.com>...
> On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> > I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
> >
> > spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
> >
> > > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > > >What do you mean by "irregular" time-spacing? If you can characterize the
> >
> > > snip
> >
> > > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > >> Hi there!
> >
> > > >> I have a signal, y which has irregular time-spacing. I would like to find
> > > >> the power spectrum of this signal. Since the signal has irregular
> > > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > > >> in Matlab?? If so, what method should I used in order to find the power
> > > >> spectrum of y?
> >
> > > snip
> >
> > > from my annotations:
> > > \begin{citation}
> > > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > > chapter will instruct you in using a periodogram-approach (called the
> > > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > > welcome to download it and give it a whirl.
> > > Regards,
> > > Brett Shoelson
> > > \end{citation}
> >
> > > \begin{citation}
> > > Luis <lfar...@uvigo.es> wrote:
> > > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > > If so, is it the same as the normal Fourier Transform or do I have to
> > > > account for anything more?
> >
> > > There are some articles published by A. Dutt an V. Rokhlin which deal
> > > with this problem. For a start you can try
> >
> > > \bibitem{dutt1}
> > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
> >
> > > \bibitem{dutt2}
> > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > > 85--100.
> >
> > > These articles contain algorithms for different types of FTs with non-
> > > equispaced data. The implementation of these algorithms isn't trivial,
> > > but they do work.
> > > \end{citation}
> >
> > > \begin{citation}
> > > we are pleased to announce the availability of Version 1.0 of a
> > > C library for computing the Nonequispaced Discrete Fourier Transform
> > > (NDFT) in one, two or three dimensions.
> > > Other common names for NFFT are non-uniform fast Fourier transform
> > > (NUFFT), generalized fast Fourier transform (GFFT),
> > > unequally-spaced fast Fourier transform (USFFT),
> > > non-equispaced fast Fourier transform (NFFT), or gridding.
> >
> > > Our library is free software based on FFTW.
> >
> > > Visit our NFFT web-page
> >
> > > http://www.math.uni-luebeck.de/potts/nfft
> >
> > > for software, documentation, and related links.
> >
> > > For further information please contact
> > > Daniel Potts (po...@math.uni-luebeck.de).
> > > \end{citation}
> >
> > > hth
> > > peter
>
> You might find the more recent discussions on this here in comp.soft-
> sys.matlab in April-May of 2008 interesting. Useful is a different
> matter.
>
> Dale B. Dalrymple


Thanks!

Guj

"guj " <gulatiakshay@gmail.com> wrote in message <h20hl6$j4r$1@fred.mathworks.com>...
> dbd <dbd@ieee.org> wrote in message <9d1ee7a3-71c6-405f-87a6-0a7fddb46b68@j9g2000prh.googlegroups.com>...
> > On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> > > I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
> > >
> > > spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
> > >
> > > > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > > > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > > > >What do you mean by "irregular" time-spacing? If you can characterize the
> > >
> > > > snip
> > >
> > > > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > > > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > > >> Hi there!
> > >
> > > > >> I have a signal, y which has irregular time-spacing. I would like to find
> > > > >> the power spectrum of this signal. Since the signal has irregular
> > > > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > > > >> in Matlab?? If so, what method should I used in order to find the power
> > > > >> spectrum of y?
> > >
> > > > snip
> > >
> > > > from my annotations:
> > > > \begin{citation}
> > > > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > > > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > > > chapter will instruct you in using a periodogram-approach (called the
> > > > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > > > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > > > welcome to download it and give it a whirl.
> > > > Regards,
> > > > Brett Shoelson
> > > > \end{citation}
> > >
> > > > \begin{citation}
> > > > Luis <lfar...@uvigo.es> wrote:
> > > > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > > > If so, is it the same as the normal Fourier Transform or do I have to
> > > > > account for anything more?
> > >
> > > > There are some articles published by A. Dutt an V. Rokhlin which deal
> > > > with this problem. For a start you can try
> > >
> > > > \bibitem{dutt1}
> > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
> > >
> > > > \bibitem{dutt2}
> > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > > > 85--100.
> > >
> > > > These articles contain algorithms for different types of FTs with non-
> > > > equispaced data. The implementation of these algorithms isn't trivial,
> > > > but they do work.
> > > > \end{citation}
> > >
> > > > \begin{citation}
> > > > we are pleased to announce the availability of Version 1.0 of a
> > > > C library for computing the Nonequispaced Discrete Fourier Transform
> > > > (NDFT) in one, two or three dimensions.
> > > > Other common names for NFFT are non-uniform fast Fourier transform
> > > > (NUFFT), generalized fast Fourier transform (GFFT),
> > > > unequally-spaced fast Fourier transform (USFFT),
> > > > non-equispaced fast Fourier transform (NFFT), or gridding.
> > >
> > > > Our library is free software based on FFTW.
> > >
> > > > Visit our NFFT web-page
> > >
> > > > http://www.math.uni-luebeck.de/potts/nfft
> > >
> > > > for software, documentation, and related links.
> > >
> > > > For further information please contact
> > > > Daniel Potts (po...@math.uni-luebeck.de).
> > > > \end{citation}
> > >
> > > > hth
> > > > peter
> >
> > You might find the more recent discussions on this here in comp.soft-
> > sys.matlab in April-May of 2008 interesting. Useful is a different
> > matter.
> >
> > Dale B. Dalrymple
>
>
> Thanks!
>
> Guj



Hi!

I never seen any NFFT Algorithm working, only on papers :) i have seen so many

Guj

"guj " <gulatiakshay@gmail.com> wrote in message <h20rea$6ls$1@fred.mathworks.com>...
> "guj " <gulatiakshay@gmail.com> wrote in message <h20hl6$j4r$1@fred.mathworks.com>...
> > dbd <dbd@ieee.org> wrote in message <9d1ee7a3-71c6-405f-87a6-0a7fddb46b68@j9g2000prh.googlegroups.com>...
> > > On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> > > > I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
> > > >
> > > > spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
> > > >
> > > > > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > > > > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > > > > >What do you mean by "irregular" time-spacing? If you can characterize the
> > > >
> > > > > snip
> > > >
> > > > > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > > > > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > > > >> Hi there!
> > > >
> > > > > >> I have a signal, y which has irregular time-spacing. I would like to find
> > > > > >> the power spectrum of this signal. Since the signal has irregular
> > > > > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > > > > >> in Matlab?? If so, what method should I used in order to find the power
> > > > > >> spectrum of y?
> > > >
> > > > > snip
> > > >
> > > > > from my annotations:
> > > > > \begin{citation}
> > > > > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > > > > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > > > > chapter will instruct you in using a periodogram-approach (called the
> > > > > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > > > > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > > > > welcome to download it and give it a whirl.
> > > > > Regards,
> > > > > Brett Shoelson
> > > > > \end{citation}
> > > >
> > > > > \begin{citation}
> > > > > Luis <lfar...@uvigo.es> wrote:
> > > > > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > > > > If so, is it the same as the normal Fourier Transform or do I have to
> > > > > > account for anything more?
> > > >
> > > > > There are some articles published by A. Dutt an V. Rokhlin which deal
> > > > > with this problem. For a start you can try
> > > >
> > > > > \bibitem{dutt1}
> > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
> > > >
> > > > > \bibitem{dutt2}
> > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > > > > 85--100.
> > > >
> > > > > These articles contain algorithms for different types of FTs with non-
> > > > > equispaced data. The implementation of these algorithms isn't trivial,
> > > > > but they do work.
> > > > > \end{citation}
> > > >
> > > > > \begin{citation}
> > > > > we are pleased to announce the availability of Version 1.0 of a
> > > > > C library for computing the Nonequispaced Discrete Fourier Transform
> > > > > (NDFT) in one, two or three dimensions.
> > > > > Other common names for NFFT are non-uniform fast Fourier transform
> > > > > (NUFFT), generalized fast Fourier transform (GFFT),
> > > > > unequally-spaced fast Fourier transform (USFFT),
> > > > > non-equispaced fast Fourier transform (NFFT), or gridding.
> > > >
> > > > > Our library is free software based on FFTW.
> > > >
> > > > > Visit our NFFT web-page
> > > >
> > > > > http://www.math.uni-luebeck.de/potts/nfft
> > > >
> > > > > for software, documentation, and related links.
> > > >
> > > > > For further information please contact
> > > > > Daniel Potts (po...@math.uni-luebeck.de).
> > > > > \end{citation}
> > > >
> > > > > hth
> > > > > peter
> > >
> > > You might find the more recent discussions on this here in comp.soft-
> > > sys.matlab in April-May of 2008 interesting. Useful is a different
> > > matter.
> > >
> > > Dale B. Dalrymple
> >
> >
> > Thanks!
> >
> > Guj
>
>
>
> Hi!
>
> I never seen any NFFT Algorithm working, only on papers :) i have seen so many
>
> Guj


Here is NDFT code written by stefan kunis and Daniel Potts, for evaluation of a trignometric polynomial at non equispaced space

function b=ndft(a,x,N,ft_opt,tflag)
%
% Computes the nonequispaced dft and its adjoint.
%
% b=ndft(a,x,N,options,tflag)
%
% a Fourier coefficients or samples
% x nodes \subset [-1/2,1/2)
% N polynomial degree
% options .method 'direct' computes the columns of the nonequispaced
% Fourier matrix by direct calls of exp
% 'horner' avoids direct calls
% 'matrix' one matrix valued call of exp
%
% tflag 'notransp' evaluate trigonometric polynomial
% 'transp' adjoint transform
%


options=ft_opt.method;

M=length(x);

freq=(-N/2):(N/2-1);

if strcmp(tflag,'notransp')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nonequispaced dft
%
% N/2-1
% f(j) = sum f_hat(k+N/2+1)*exp(-2*pi*i*k*x(j)), 1 <= j <= M.
% k=-N/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  f_hat=a;
  switch options
    case 'direct'
      f=zeros(M,1);
      for k=1:N
        f=f+f_hat(k).*exp(-2*pi*i*x*freq(k));
      end;
    case 'horner'
      f=zeros(M,1);
      exp_x=exp(2*pi*i*x);
      for k=1:N
        f=(f+f_hat(k)).*exp_x;
      end;
      f=f.*exp(-N*pi*i*x);
    case 'matrix'
      f=exp(-2*pi*i*x*freq)*f_hat;
    otherwise
      disp('unknown option')
  end;
  b=f;

elseif strcmp(tflag,'transp')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% adjoint nonequispaced dft
%
% M
% f_hat(k+N/2+1) = sum f(j)*exp(2*pi*i*k*x(j)), -N/2 <= k <= N/2-1.
% j=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  f=a;
  switch options
    case 'direct'
      f_hat=zeros(N,1);
      for k=1:N
        f_hat(k)=exp(-2*pi*i*x*freq(k))' *f;
        % exp(-2*pi*i*freq(k)*x)' *f differs by 10^-11!?
      end;
    case 'horner'
      f_hat=zeros(N,1);
      f=f.*exp(-N*pi*i*x);
      exp_x=exp(2*pi*i*x);
      for k=1:N
        f_hat(k)=sum(f);
        f=f.*exp_x;
      end;
    case 'matrix'
      f_hat=(exp(-2*pi*i*x*freq)') *f;
    otherwise
      disp('unknown option')
  end;
  b=f_hat;

end;



------------------------------------------------------------------------------------
1. Whats the reason to put it exactly between -1/2 to 1/2
2. And, I never understand the mathematicians why they started with frequency to time rather than time to frequency any special reasons for that
3. Plz check this code for uniform sampling with conventional FFT..According to theory it should give the same result, not for me

"guj " <gulatiakshay@gmail.com> wrote in message <h3frp4$5m3$1@fred.mathworks.com>...
> "guj " <gulatiakshay@gmail.com> wrote in message <h20rea$6ls$1@fred.mathworks.com>...
> > "guj " <gulatiakshay@gmail.com> wrote in message <h20hl6$j4r$1@fred.mathworks.com>...
> > > dbd <dbd@ieee.org> wrote in message <9d1ee7a3-71c6-405f-87a6-0a7fddb46b68@j9g2000prh.googlegroups.com>...
> > > > On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> > > > > I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
> > > > >
> > > > > spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
> > > > >
> > > > > > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > > > > > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > > > > > >What do you mean by "irregular" time-spacing? If you can characterize the
> > > > >
> > > > > > snip
> > > > >
> > > > > > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > > > > > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > > > > >> Hi there!
> > > > >
> > > > > > >> I have a signal, y which has irregular time-spacing. I would like to find
> > > > > > >> the power spectrum of this signal. Since the signal has irregular
> > > > > > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > > > > > >> in Matlab?? If so, what method should I used in order to find the power
> > > > > > >> spectrum of y?
> > > > >
> > > > > > snip
> > > > >
> > > > > > from my annotations:
> > > > > > \begin{citation}
> > > > > > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > > > > > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > > > > > chapter will instruct you in using a periodogram-approach (called the
> > > > > > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > > > > > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > > > > > welcome to download it and give it a whirl.
> > > > > > Regards,
> > > > > > Brett Shoelson
> > > > > > \end{citation}
> > > > >
> > > > > > \begin{citation}
> > > > > > Luis <lfar...@uvigo.es> wrote:
> > > > > > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > > > > > If so, is it the same as the normal Fourier Transform or do I have to
> > > > > > > account for anything more?
> > > > >
> > > > > > There are some articles published by A. Dutt an V. Rokhlin which deal
> > > > > > with this problem. For a start you can try
> > > > >
> > > > > > \bibitem{dutt1}
> > > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
> > > > >
> > > > > > \bibitem{dutt2}
> > > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > > > > > 85--100.
> > > > >
> > > > > > These articles contain algorithms for different types of FTs with non-
> > > > > > equispaced data. The implementation of these algorithms isn't trivial,
> > > > > > but they do work.
> > > > > > \end{citation}
> > > > >
> > > > > > \begin{citation}
> > > > > > we are pleased to announce the availability of Version 1.0 of a
> > > > > > C library for computing the Nonequispaced Discrete Fourier Transform
> > > > > > (NDFT) in one, two or three dimensions.
> > > > > > Other common names for NFFT are non-uniform fast Fourier transform
> > > > > > (NUFFT), generalized fast Fourier transform (GFFT),
> > > > > > unequally-spaced fast Fourier transform (USFFT),
> > > > > > non-equispaced fast Fourier transform (NFFT), or gridding.
> > > > >
> > > > > > Our library is free software based on FFTW.
> > > > >
> > > > > > Visit our NFFT web-page
> > > > >
> > > > > > http://www.math.uni-luebeck.de/potts/nfft
> > > > >
> > > > > > for software, documentation, and related links.
> > > > >
> > > > > > For further information please contact
> > > > > > Daniel Potts (po...@math.uni-luebeck.de).
> > > > > > \end{citation}
> > > > >
> > > > > > hth
> > > > > > peter
> > > >
> > > > You might find the more recent discussions on this here in comp.soft-
> > > > sys.matlab in April-May of 2008 interesting. Useful is a different
> > > > matter.
> > > >
> > > > Dale B. Dalrymple
> > >
> > >
> > > Thanks!
> > >
> > > Guj
> >
> >
> >
> > Hi!
> >
> > I never seen any NFFT Algorithm working, only on papers :) i have seen so many
> >
> > Guj
>
>
> Here is NDFT code written by stefan kunis and Daniel Potts, for evaluation of a trignometric polynomial at non equispaced space
>
> function b=ndft(a,x,N,ft_opt,tflag)
> %
> % Computes the nonequispaced dft and its adjoint.
> %
> % b=ndft(a,x,N,options,tflag)
> %
> % a Fourier coefficients or samples
> % x nodes \subset [-1/2,1/2)
> % N polynomial degree
> % options .method 'direct' computes the columns of the nonequispaced
> % Fourier matrix by direct calls of exp
> % 'horner' avoids direct calls
> % 'matrix' one matrix valued call of exp
> %
> % tflag 'notransp' evaluate trigonometric polynomial
> % 'transp' adjoint transform
> %
>
>
> options=ft_opt.method;
>
> M=length(x);
>
> freq=(-N/2):(N/2-1);
>
> if strcmp(tflag,'notransp')
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % nonequispaced dft
> %
> % N/2-1
> % f(j) = sum f_hat(k+N/2+1)*exp(-2*pi*i*k*x(j)), 1 <= j <= M.
> % k=-N/2
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> f_hat=a;
> switch options
> case 'direct'
> f=zeros(M,1);
> for k=1:N
> f=f+f_hat(k).*exp(-2*pi*i*x*freq(k));
> end;
> case 'horner'
> f=zeros(M,1);
> exp_x=exp(2*pi*i*x);
> for k=1:N
> f=(f+f_hat(k)).*exp_x;
> end;
> f=f.*exp(-N*pi*i*x);
> case 'matrix'
> f=exp(-2*pi*i*x*freq)*f_hat;
> otherwise
> disp('unknown option')
> end;
> b=f;
>
> elseif strcmp(tflag,'transp')
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % adjoint nonequispaced dft
> %
> % M
> % f_hat(k+N/2+1) = sum f(j)*exp(2*pi*i*k*x(j)), -N/2 <= k <= N/2-1.
> % j=1
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> f=a;
> switch options
> case 'direct'
> f_hat=zeros(N,1);
> for k=1:N
> f_hat(k)=exp(-2*pi*i*x*freq(k))' *f;
> % exp(-2*pi*i*freq(k)*x)' *f differs by 10^-11!?
> end;
> case 'horner'
> f_hat=zeros(N,1);
> f=f.*exp(-N*pi*i*x);
> exp_x=exp(2*pi*i*x);
> for k=1:N
> f_hat(k)=sum(f);
> f=f.*exp_x;
> end;
> case 'matrix'
> f_hat=(exp(-2*pi*i*x*freq)') *f;
> otherwise
> disp('unknown option')
> end;
> b=f_hat;
>
> end;
>
>
>
> ------------------------------------------------------------------------------------
> 1. Whats the reason to put it exactly between -1/2 to 1/2
    1. Well, i think its just because its periodic and in this paper they used this for truncating their gaussian later so thats the number they have chosen
> 2. And, I never understand the mathematicians why they started with frequency to time rather than time to frequency any special reasons for that
> 3. Plz check this code for uniform sampling with conventional FFT..According to theory it should give the same result, not for me






Rapid Computation of the Discrete Fourier Transform by chirs anderson 1996 ...SIAM Journal of scintific computing...can any one mail me this paper at pinkfloydindia@yahoo.com
I have only access to paper after 1997 in my university

"guj " <gulatiakshay@gmail.com> wrote in message <h3g81f$n3k$1@fred.mathworks.com>...
> "guj " <gulatiakshay@gmail.com> wrote in message <h3frp4$5m3$1@fred.mathworks.com>...
> > "guj " <gulatiakshay@gmail.com> wrote in message <h20rea$6ls$1@fred.mathworks.com>...
> > > "guj " <gulatiakshay@gmail.com> wrote in message <h20hl6$j4r$1@fred.mathworks.com>...
> > > > dbd <dbd@ieee.org> wrote in message <9d1ee7a3-71c6-405f-87a6-0a7fddb46b68@j9g2000prh.googlegroups.com>...
> > > > > On Jun 25, 7:53 am, "Akshay " <gulatiaks...@gmail.com> wrote:
> > > > > > I dont think, the matlab code given on ur site works. It gives crappy results. Still have to try the c library. Matlab code should give results, in c its all about performance. Right now i am only concern whether algorithm works or not. It seems impossible to reconstruct the original data.
> > > > > >
> > > > > > spellu...@mathematik.tu-darmstadt.de (Peter Spellucci) wrote in message <apr8nf$q0...@fb04373.mathematik.tu-darmstadt.de>...
> > > > > >
> > > > > > > In article <1R3w9.38629$iV1.8...@nwrddc02.gnilink.net>,
> > > > > > > "Michael Pender" <michael.pen...@nanochron.com> writes:
> > > > > > > >What do you mean by "irregular" time-spacing? If you can characterize the
> > > > > >
> > > > > > > snip
> > > > > >
> > > > > > > >Chia C Chong <Chia.Ch...@ee.ed.ac.uk> wrote in message
> > > > > > > >news:ap102a$6ce$1@scotsman.ed.ac.uk...
> > > > > > > >> Hi there!
> > > > > >
> > > > > > > >> I have a signal, y which has irregular time-spacing. I would like to find
> > > > > > > >> the power spectrum of this signal. Since the signal has irregular
> > > > > > > >> time-spacing, I guess I can't simply perform the DFT using the fft command
> > > > > > > >> in Matlab?? If so, what method should I used in order to find the power
> > > > > > > >> spectrum of y?
> > > > > >
> > > > > > > snip
> > > > > >
> > > > > > > from my annotations:
> > > > > > > \begin{citation}
> > > > > > > You might take a look at "Spectral Analysis of Unevenly Sampled Data"
> > > > > > > in any of the Numerical Recipes books (Press, Teukolsky, etc.). The
> > > > > > > chapter will instruct you in using a periodogram-approach (called the
> > > > > > > Lomb-Scargle Periodogram) to analyze your unevenly sampled data. I've
> > > > > > > posted a MATLAB version of Lomb-Scargle on MATLAB Central; you're
> > > > > > > welcome to download it and give it a whirl.
> > > > > > > Regards,
> > > > > > > Brett Shoelson
> > > > > > > \end{citation}
> > > > > >
> > > > > > > \begin{citation}
> > > > > > > Luis <lfar...@uvigo.es> wrote:
> > > > > > > > Can I do a Fourier Transform of data thats not equaly spaced in time?
> > > > > > > > If so, is it the same as the normal Fourier Transform or do I have to
> > > > > > > > account for anything more?
> > > > > >
> > > > > > > There are some articles published by A. Dutt an V. Rokhlin which deal
> > > > > > > with this problem. For a start you can try
> > > > > >
> > > > > > > \bibitem{dutt1}
> > > > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > > > Data}, SIAM J. Sci. Comp. {\bf 14} (1993), 1368--1393.
> > > > > >
> > > > > > > \bibitem{dutt2}
> > > > > > > A.~Dutt, V.~Rokhlin, {\it Fast Fourier Transforms for Nonequispaced
> > > > > > > Data, II}, Applied and Computational Harmonic Analysis {\bf 2} (1995),
> > > > > > > 85--100.
> > > > > >
> > > > > > > These articles contain algorithms for different types of FTs with non-
> > > > > > > equispaced data. The implementation of these algorithms isn't trivial,
> > > > > > > but they do work.
> > > > > > > \end{citation}
> > > > > >
> > > > > > > \begin{citation}
> > > > > > > we are pleased to announce the availability of Version 1.0 of a
> > > > > > > C library for computing the Nonequispaced Discrete Fourier Transform
> > > > > > > (NDFT) in one, two or three dimensions.
> > > > > > > Other common names for NFFT are non-uniform fast Fourier transform
> > > > > > > (NUFFT), generalized fast Fourier transform (GFFT),
> > > > > > > unequally-spaced fast Fourier transform (USFFT),
> > > > > > > non-equispaced fast Fourier transform (NFFT), or gridding.
> > > > > >
> > > > > > > Our library is free software based on FFTW.
> > > > > >
> > > > > > > Visit our NFFT web-page
> > > > > >
> > > > > > > http://www.math.uni-luebeck.de/potts/nfft
> > > > > >
> > > > > > > for software, documentation, and related links.
> > > > > >
> > > > > > > For further information please contact
> > > > > > > Daniel Potts (po...@math.uni-luebeck.de).
> > > > > > > \end{citation}
> > > > > >
> > > > > > > hth
> > > > > > > peter
> > > > >
> > > > > You might find the more recent discussions on this here in comp.soft-
> > > > > sys.matlab in April-May of 2008 interesting. Useful is a different
> > > > > matter.
> > > > >
> > > > > Dale B. Dalrymple
> > > >
> > > >
> > > > Thanks!
> > > >
> > > > Guj
> > >
> > >
> > >
> > > Hi!
> > >
> > > I never seen any NFFT Algorithm working, only on papers :) i have seen so many
> > >
> > > Guj
> >
> >
> > Here is NDFT code written by stefan kunis and Daniel Potts, for evaluation of a trignometric polynomial at non equispaced space
> >
> > function b=ndft(a,x,N,ft_opt,tflag)
> > %
> > % Computes the nonequispaced dft and its adjoint.
> > %
> > % b=ndft(a,x,N,options,tflag)
> > %
> > % a Fourier coefficients or samples
> > % x nodes \subset [-1/2,1/2)
> > % N polynomial degree
> > % options .method 'direct' computes the columns of the nonequispaced
> > % Fourier matrix by direct calls of exp
> > % 'horner' avoids direct calls
> > % 'matrix' one matrix valued call of exp
> > %
> > % tflag 'notransp' evaluate trigonometric polynomial
> > % 'transp' adjoint transform
> > %
> >
> >
> > options=ft_opt.method;
> >
> > M=length(x);
> >
> > freq=(-N/2):(N/2-1);
> >
> > if strcmp(tflag,'notransp')
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> > % nonequispaced dft
> > %
> > % N/2-1
> > % f(j) = sum f_hat(k+N/2+1)*exp(-2*pi*i*k*x(j)), 1 <= j <= M.
> > % k=-N/2
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> >
> > f_hat=a;
> > switch options
> > case 'direct'
> > f=zeros(M,1);
> > for k=1:N
> > f=f+f_hat(k).*exp(-2*pi*i*x*freq(k));
> > end;
> > case 'horner'
> > f=zeros(M,1);
> > exp_x=exp(2*pi*i*x);
> > for k=1:N
> > f=(f+f_hat(k)).*exp_x;
> > end;
> > f=f.*exp(-N*pi*i*x);
> > case 'matrix'
> > f=exp(-2*pi*i*x*freq)*f_hat;
> > otherwise
> > disp('unknown option')
> > end;
> > b=f;
> >
> > elseif strcmp(tflag,'transp')
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> > % adjoint nonequispaced dft
> > %
> > % M
> > % f_hat(k+N/2+1) = sum f(j)*exp(2*pi*i*k*x(j)), -N/2 <= k <= N/2-1.
> > % j=1
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> >
> > f=a;
> > switch options
> > case 'direct'
> > f_hat=zeros(N,1);
> > for k=1:N
> > f_hat(k)=exp(-2*pi*i*x*freq(k))' *f;
> > % exp(-2*pi*i*freq(k)*x)' *f differs by 10^-11!?
> > end;
> > case 'horner'
> > f_hat=zeros(N,1);
> > f=f.*exp(-N*pi*i*x);
> > exp_x=exp(2*pi*i*x);
> > for k=1:N
> > f_hat(k)=sum(f);
> > f=f.*exp_x;
> > end;
> > case 'matrix'
> > f_hat=(exp(-2*pi*i*x*freq)') *f;
> > otherwise
> > disp('unknown option')
> > end;
> > b=f_hat;
> >
> > end;
> >
> >
> >
> > ------------------------------------------------------------------------------------
> > 1. Whats the reason to put it exactly between -1/2 to 1/2
> 1. Well, i think its just because its periodic and in this paper they used this for truncating their gaussian later so thats the number they have chosen
> > 2. And, I never understand the mathematicians why they started with frequency to time rather than time to frequency any special reasons for that
> > 3. Plz check this code for uniform sampling with conventional FFT..According to theory it should give the same result, not for me
>
>
>
>
>
>
> Rapid Computation of the Discrete Fourier Transform by chirs anderson 1996 ...SIAM Journal of scintific computing...can any one mail me this paper at pinkfloydindia@yahoo.com
> I have only access to paper after 1997 in my university









Hello All! I found the tutorial of NDFT

http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NUFT/node7.html#SECTION00032000000000000000


Although I am not able to reproduce the result as they had done for NDFT here, if any body gets the same result for NDFT please let me know, do try it waiting for your comments

If you are still looking for code, I have submitted some 1D/2D/3D NUFFT Matlab functions that implement "Fast Gaussian Gridding."