Dft for Nonuniform Spacing
Greg Heath
Comments and corrections are welcome.
This topic has been discussed many times over the
years in comp.soft-sys.matlab, comp.dsp, and
sci.math.num-analysis. However, I have not found
a good basic online explanation in a single source.
For example in
Newsgroups: comp.dsp, comp.soft-sys.matlab
Date: 23 Oct 2002 23:27:21 -0700
Subject: Re: DFT of an irregular time-spacing signal
http://groups.google.com/group/...
comp.soft-sys.matlab/msg/39f42567bbd8374b
I proposed a nonuniform sampling generalization of
AJ Johnson's MATLAB dft function. Several
questions were raised (but not satisfactorily
answered!) regarding sufficient conditions on
max(f) and min(df)for an accurate reconstruction
of the original nonuniformly sampled function x(t).
First a clarification of terms
DFT Any finite difference approximation to the
finite-domain Fourier integral
X(f) = INT(t=0,T){dt exp(-2*pi*j*t) x(t)}
FFT Special, fast, implementation of the dft
for a N point uniform sampling of x to be
evaluated at N uniformly spaced frequencies.
In a recent post
Newsgroups: comp.soft-sys.matlab
Date: Thu, 17 Apr 2008 14:28:31 -0700 (PDT)
Subject: Re: Spectra of Unequally Sampled Data
http://groups.google.com/group/...
comp.soft-sys.matlab/msg/63159e7825b4e96e
I provided details by using a piecewise constant
(PWC) approximation of x and represented the
integral as a sum of rectangular areas. A more
accurate implentation is easily obtained by using
a piecewise linear (PWL) approximation of x and
representing the integral as a sum of trapezoidal
areas. Although higher order interpolation can
be used, it won't be considered in this discussion.
In the following, t, x, f and X are assumed to be
column vectors. The results of the PWC and PWL
aprroximations are
N = length(t);
dt = diff(t);
dts = ([0 dt]+[dt 0])/2;
W = exp(-2*pi*j* f*t');
X = W*(x.*dts);
and
tm = t(1:end-1)+t(2:end))/2;
xm = x(1:end-1)+x(2:end))/2;
Wm = exp(-2*pi*j* f*tm');
X = Wm*(xm.*dt);
respectively.
Note that the lengths of tm and xm are N-1 and
x can be recovered from xm only if sufficient
apriori information about the original points
is known. For example, if the mean of x is zero,
x(1) = x1; % x1 is arbitrary
for i = 2:N
x(i) = 2*xm(i-1)-x(i-1);
end
x = x - mean(x);
The corresponding Fourier inversion formulae are
(1) x = W'*(X.*dfs); (2) xm = Wm'*(X.*dfs);
( For the purpose of clarity, fm and Xm will not
be considered here). However, (1) and (2) are not
valid when either t or f are nonuniformly spaced
because W and Wm are no longer unitary (e.g.,
W'*W ~= I).
Although the least square approximations
(3) x = (W\X)./dts; (4) xm = (Wm\X)./dt;
could be used, it is more satisfying to start with
the explicit summation over sinusoids in (1),(2),
and derive the least square approximations
(5) X = (W'\x)./dfs; (6) X = (Wm'\xm)./dfs;
If min(dt) is too small, two or more columns of
W will be highly correlated and W will be
ill-conditioned. Similarly, if min(df) is too
small, W will be ill-conditioned because two or
more rows of W will be highly correlated.
Therefore it is reasonable to assume that
min(df)* min(dt) >= C = constant
is a necessary condition for accurate
representations of x and it's spectrum. Although
df*dt = 1/N for uniform sampling when M =
length(f) = length(x) = N, it is not obvious what
C should be when sampling is nonuniform and/or
the lengths of f and t are not equal. However,
the outer product constraint
1/sqrt(M*N) <= df*dt'
seems to be reasonable.
Consider the trivial case
t = [t1 t2], dt > 0,
f = [f1 f2], df > 0,
X(f)= x1*exp(-2*pi*j*t1*f)+ x2*exp(-2*pi*j*t2*f).
It is trivial to show that x1 and x2 can be
estimated from samples at f if
exp(2*pi*j*df*dt) ~= 0
i.e.,
df * dt ~= 0,1,2,...
Therefore, it is reasonable to consider
the outer product matrix bounds
1/sqrt(M*N) <= df * dt' < 1
as sufficient.
For the familiar case of uniform sampling x and
X are periodic with M = N, df*dt = 1/N, and
df = 1/((max(t)+dt),
max(f) = (1/dt)-df = [max(t)/(max(t)+dt)]*(1/dt).
The case of nonuniform sampling can be interpreted
using the following model of uniform sampling:
Define N0, dt0 and the integer sequence
n = [n1,n2,...,nN]' so that that dt0 is the largest
interval for which
t(i) = n(i)*dt0
max(t) = (N0-1)*dt0.
Now consider resampling x at the rate 1/dt0 to obtain
t0(j) = (j-1)*dt0, j = 1,2,...N0
for i = 1:N
if t0(j) = t(i)
x0(j) = x(i);
else
x0(j) = 0;
end
end.
Then x0 is identical to x at sampling points and zero
in between.
Since (N0-1)*dt0 = max(t0) = max(t)= (N-1)*dt, the
corresponding transform, X0, is defined on f0 with
df0 = 1/(max(t0)+dt0)
= ((N0-1)/(N-1))*(N/N0)*df
~ df
max(f0) = [max(t0)/(max(t0)+dt0)]*(1/dt0).
= ((N0-1)/(N-1))*max(f)
~ (N0/N)*max(f)
Typically, df0 is only slightly smaller than df
but max(f0) can be orders of magnitude greater
than max(f)!
The above result is not magic. Nonuniform sampling
sampling can significantly modify the periodic
property that allows high frequency information to
leak into lower frequencies though the phenomenon
of aliasing.
A specific nonuniform sampling pattern dt defines
N0 but df and max(f) are at the discretion of the
user.
The bottom line is that it is not unreasonable to
make the choices df and M so that
df = 1/( max(t)+ min(dt))
(M-1)*df = max(f) <= 1/min(dt) < M*df
Hope this helps.
Greg
dbd
Helps with what Greg?
The question that triggered the posts you reference from 2002 was:
> 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.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Is there an answer yet? If so, what is it?
Dale B. Dalrymple
http://dbdimages.com
DigitalSignal
detail, could you let us know what is the conclusion? Did you find a
fast transform algorithm to map the frequency-varying signal into
spectrum domain?
James
www.go-ci.com
Greg Heath
Sorry for the misdirection. I was only concerned with
three questions.
Q1. What formula most accurately reconstructs the
original signal?
A1. a. If the spectrum is obtained from the least squares
formula (5) X = (W'\x)./dfs, then the reconstruction
(1) x = W'*X.*dfs will be more accurate than when the
spectrum is obtained from the usual Fourier formula
X = W*(x.*dts).
b. However, if the spectrum is obtained from the usual
Fourier formula X = W*(x.*dts), then the least squares
reconstruction formula (3) x = (W\X)./dts will be more
accurate than the Fourier formula (1) x = W'*X.*dfs.
Q2. What is a reasonable choice for df?
A2. df = 1/(max(t)+ min(dt)) is sufficient.
Q3. If f = df*(0:M-1), what is a reasonable choice for M?
A3. M-1 <= 1/min(df*dt) < M
In response to the original 2002 question, the above the LS
expressions
represent linear equations of the form A*u = v, with minimum square
error solutions u = B*v on which perturbation analyses have been done
in the literature provided A is of full rank. When A is rank
deficient, unique
solutions probably have to be guaranteed imposing additional
constraints.
This is standard LS theory. However, I don't have any references handy
...Sorry.
Hope this helps,
Greg
Greg Heath
> I am very interested in this topic. Before I read the program in
> detail, could you let us know what is the conclusion?
See my previous reply.
> Did you find a fast transform algorithm to map the frequency-varying
> signal into spectrum domain?
Almost. In the second post that I referenced I compared my dft with
MATLABs
fft for the uniformly sampled function exp(-t/T).*cos(t), 0 <= t <= T
= 4*pi
and found that the fft is only 500 to 1000 times faster when N = 128.
Therefore
I only have to do a little more tweaking.
Seriously, there are fast transforms for nonuniformly sampled data. I
am
not familiar with them. The search keyword is NFFT...and there may be
free downloads available in MATLAB and/or C.
Obviously, my curiosity of this bandied about topic has just scratched
the
surface.
Hope this helps.
Greg
dbd
> ....
> ...
The question was and still is:
> How closely does an evaluation of the spectrum information reproduce the
> original data?
> Jerry
> --
What numerical evaluation and comparison have you made to demonstrate
the validity of your subjective claims of 'reasonableness' against the
other reconstructions? What level of deviation have you used for your
decisions of 'reasonable'?
Q1 is to the point but unevaluated. Is the difference significant and
does it matter?.
Q2 and Q3 include on undefined terminology. Numbers in A1 would let
the reader evaluate A2 and A3
Hoping for help,
Dale B. Dalrymple
Dieter Britz
Check Britz & Antonia, Meas. Sci. Technol. 7 (1996) 1042 and references
therein.
--
Dieter Britz (britz<at>chem.au.dk)
Greg Heath
The reason I looked into this was to satify my own curiosity
about the answer to the three basic questions stated above.
About 20 years ago a colleague who was analyzing nonuniformly
spaced Doppler radar data fron a precessing and spinning
target asked me to help him with his analysis. I know I
considered Lomb-Scargle but I honestly don't remember what I
finally recommended.
I posted this to
(1) Suggest what appears to be an improvement to the modified
AJ Johnson dft code by using trapezoidal instead of rectangular
area approximations to the Fourier integrals.
(2) Suggest LS for the inverse time <--> frequency transformation.
(3) Give reasonable criteria for choosing the sampling frequencies.
(4) See if I was on the right track by inviting comments from
those with more experience.
#2 was suggested by several others in old threads. I could not find
a quick answer to #3 via a Google search. As for #1, of course higher
order interpolations are available. However, they seem to be frowned
on by experts because of the high frequency artifacts they can
introduce. Furthermore, using higher order interpolation was not
going to help me answer my 3 basic questions.
As for #4, programming experiments would be a waste of time if my
reasoning were faulty.
> What numerical evaluation and comparison have you made to demonstrate
> the validity of your subjective claims of 'reasonableness' against the
> other reconstructions?
None. See my last comment.
> What level of deviation have you used for your
> decisions of 'reasonable'?
"Reasonable" is subjective. Unless there is an obvious flaw in my
reasoning, parameterized computer experiments would help quantify my
qualitative conclusions.
> Q1 is to the point but unevaluated. Is the difference significant and
> does it matter?.
It all depends on the underlying function and the sampling pattern.
The points are (1) it's obviously a better approximation to the
integral with a minor change in code.(2) a low order interpolation
is less likely to cause high frequency artifacts.
> Q2 and Q3 include on undefined terminology.
Not sure what that means.
> Numbers in A1 would let
> the reader evaluate A2 and A3
Agree.
> Hoping for help,
>
> Dale B. Dalrymple
Thanks, again, for your previous help. I had been recommending
AJJs dft code to others who needed code for nonuniform sampling.
Greg
dbd
When you play this game of numerical approximation I don't know how
you can discuss 'reasonableness' without numerical evaluation. It
serves two purposes. 1) It forces a selection of a 'reasonable'
instantiation of nonuniform sampling. This is something you have
failed to do which prevents any rational response to your ponderings.
How can we decide what is a reasonable input? 2) It evaluates how well
the processing achieves it's goal by presenting error values that
allow us to apply our (probably different) interpretations of
reasonable error to a known referent.
Every day people apply uniform algorithms to data collected with
imperfect clocks. In many cases, perhaps most, it is a reasonable
thing to do. But manufacturers of high speed data converters have
provided applications literature to enable users to subjectively
evaluate what amount of jitter will produce a reasonable degradation
in performance.
This topic is a deliberate foray from the beaten track. Tools like
Matlab let us make quick evaluation of how far into the wilds we have
wandered. They make the algorithm description executable. Perhaps the
number of lines of Matlab code and the timing comparisons in your
postings have given us unreasonable expectations of useful algorithm
evaluation.
Dale B. Dalrymple
Greg Heath
Not sure what you mean here. Although some have purposely
used nonuniform sampling to avoid aliasing, typically
nonuniform sampling is something the user is presented
with. The most common are
(1) Piecewise uniform sampling
(2) Uniform sampling with deterministic or random gaps
of missing data.
(3) Almost periodic sampling with times "jittered" about
the uniform sampling times.
(4) Nonuniform sampling with monotonically increasing
and/or decreasing sampling times.
However, it would be useful to find out what examples
have been demonstrated in online sources.
> This is something you have failed to do which prevents
> any rational response to your ponderings.
> How can we decide what is a reasonable input?
Not sure what you mean here either. What I listed above
or input to the DFT i.e., df and M?
> 2) It evaluates how well
> the processing achieves it's goal by presenting error values that
> allow us to apply our (probably different) interpretations of
> reasonable error to a known referent.
There are too many variations of interesting functions and
sample patterns to choose from. Do you have one that you
would like to see?
> Every day people apply uniform algorithms to data collected with
> imperfect clocks. In many cases, perhaps most, it is a reasonable
> thing to do. But manufacturers of high speed data converters have
> provided applications literature to enable users to subjectively
> evaluate what amount of jitter will produce a reasonable degradation
> in performance.
If you are interested in a jitter example. I'll cook one up
for a demo. Although it shouldn't take long to do, I won't
be able to work on it for a week or so.
> This topic is a deliberate foray from the beaten track.
It's a deliberate response to those who are still asking
newsgroup questions about how to do an FFT on nonuniformly
sampled data. In particular, how to easily modify and use
AJJs MATLAB algorithm.
> Tools like
> Matlab let us make quick evaluation of how far into the wilds we have
> wandered. They make the algorithm description executable. Perhaps the
> number of lines of Matlab code and the timing comparisons in your
> postings have given us unreasonable expectations of useful algorithm
> evaluation.
Not sure what that means. My timings were on uniform sampling.
Showing that the DFT is 1000 times slower than a FFT shouldn't
raise anyones expectations.
Greg
P.S. The original question in the 2002 post was
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?
Greg Heath
> On Apr 29, 3:59 pm, dbd <d...@ieee.org> wrote:
> > On Apr 29, 4:04 am, Greg Heath <he...@alumni.brown.edu> wrote:
> > > On Apr 29, 12:27 am, dbd <d...@ieee.org> wrote:
> > > > On Apr 28, 5:05 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > This topic is a deliberate foray from the beaten track.
>
> It's a deliberate response to those who are still asking
> sampled data. In particular, how to easily modify and use
> AJJs MATLAB algorithm.
In particular the goal of this post was an explanation of how to
modify
and use AJJs MATLAB algorithm.
> > Tools like
> > Matlab let us make quick evaluation of how far into the wilds we have
> > wandered. They make the algorithm description executable. Perhaps the
> > number of lines of Matlab code and the timing comparisons in your
> > postings have given us unreasonable expectations of useful algorithm
> > evaluation.
>
> Not sure what that means. My timings were on uniform sampling.
> Showing that the DFT is 1000 times slower than aFFTshouldn't
> raise anyones expectations.
>
> Greg
>
> P.S. The original question in the 2002 post was
>
> 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 thefft
> command in Matlab?? If so, what method should I used in order to find
> the power spectrum of y?
and, the 2008 post
I need to calculate the spectra of unequally sampled data.
The data is from LDV (Laser Doppler Velocimetry). Two
packages from the Mathworks User File Exchange come to
mind, Automatic Spectra Analysis and the ARMASA toolbox.
The documentation for these is a bit tough for me to
understand.
So, given a set of data points sampled at unequal time
intervals, how do I get to spectral results?
Greg
dbd
Greg
My observation about the topic of this discussion as 'off the beaten
track' is about the nature of the topic, not the content of the
discussion. The topic keeps alive because people are interested in it.
They are interested because they want to learn about something that
they know they don't understand. I made the remark to emphasize that
the topic is one where few people have a good feeling for what is
'reasonable'. In light of that, how can you expect people in comp or
sci groups to look for anything but numbers? That's what we use to
decide what is 'reasonable'. This discussion began with a 200 line
post full of code describing your concepts. Perhaps that let me
mislead myself about how close you seemed to be to being able to give
numbers. Your efforts are appreciated and I hope to see some results.
Dale B. Dalrymple
Greg Heath
> On May 1, 4:09 am,Greg Heath<he...@alumni.brown.edu> wrote:
>
>
>
>
>
> > On May 1, 6:37 am,Greg Heath<he...@alumni.brown.edu> wrote:
>
> > > On Apr 29, 3:59 pm, dbd <d...@ieee.org> wrote:
> > > > > On Apr 29, 12:27 am, dbd <d...@ieee.org> wrote:
My goal is basic. If you use the simple algorithms discussed above,
what constraints are on fmax and constant df ? Why? What happens if
they are violated?
Fast and/or accurate algorithms already exist (e.g., search on
NFFT) .
However, the sophistication obscures the basics. .
Greg
Rathna....@gmail.com
non-uniformed spaced signal with constant magnitude(same as 1). So
basically the time sampling is non-uniform.
I would like to do fourier transforms on this time signals to see the
spectrum. I want to measure the asymmetry in these time signals. Can
you recommend me any simple method to do this?
/Best Regards/RC
Greg Heath
If you are adventurous, you can try my unguaranteed modification
of AJJs DFT algorithm included below.
Alternatively, you can download the LOMB-SCARGLE algorithm from
the MATHWORKS website or a fast DFT from the NFFT website.
Hope this helps.
Greg
function [XFT,XLS,NMSEFT,NMSELS] = DFTgh1(x,t,f)
% function [XFT,XLS,NMSEFT,NMSELS] = DFTgh1(x,t,f)
%
% Modification of AJ Johnson's dft for nonuniform sampling
%
% Computes XFT (Discrete Fourier Transform) at frequencies
% given in f, given samples x taken at times t:
%
% XFT(f) = sum(k=1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
% = W *(x.*dts)
%
% where dts is a symmetrized modification of diff(t).
%
% Also computes the Least-Squared-Error Spectrum at
% frequencies given in f, given samples x taken at
% times t:
%
% XLS(f) = (W'\x)./dfs;
%
% where dfs is a symmetrized modification of diff(f).
%
% NMSEFT is the normalized mean-square-error of reconstucting
% x from X using the Inverse Fourier Transform formula. If
% mean(x) = 0, then the MSE is unnormalized.
%
% NMSELS is the normalized mean-square-error of reconstucting
% x from X using Least Squares. If mean(x) = 0, then the MSE
% is unnormalized.
%
% For comparison with MATLAB's FFT when the spacing is uniform,
% double the end values x(1) and x(end) and divide X by dt0 =
% mean(diff(t))
x = x(:); % Format 'x' into a column vector
t = t(:); % Format 't' into a column vector
f = f(:); % Format 'f' into a column vector
N = length(x);
if length(t)~= N
error('x and t do not have the same length')
end;
dt = diff(t); % asymmetric "dt"
dts = 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
meanx = sum(x.*dts)/sum(dts);
df = diff(f); % asymmetric "df"
dfs = 0.5*([df; 0]+[0; df]); % symmetric "df"
W = exp(-2*pi*j * f*t');
XFT = W * (x.*dts);
XLS = (W'\x)./dfs;
xFT = real(W'*(XFT.*dfs));
xLS = real(W'*(XLS.*dfs));
MSE0 = mse(abs(x-meanx));
if MSE0 == 0, MSE0 = 1, end;
NMSEFT = mse(abs(x-xFT))/MSE0;
NMSELS = mse(abs(x-xLS))/MSE0;
return
Bruce Scott TOK
One alternative would be to interpolate first and then FFT. You would
have to analyse the effect of your interpolation scheme, but the
additional scales added by something like 3rd order Lagrange
interpolation should lie outside your actual resolved spectrum.
--
ciao,
Bruce
drift wave turbulence: http://www.rzg.mpg.de/~bds/