Dft for Nonuniform Spacing 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
Dft for Nonuniform Spacing
Started by ●April 28, 2008
Reply by ●April 28, 20082008-04-28
On Apr 28, 12:33 pm, Greg Heath <he...@alumni.brown.edu> wrote:> Dft for Nonuniform Spacing > > 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. > > GregHelps 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
Reply by ●April 28, 20082008-04-28
I am very interested in this topic. Before I read the program in 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
Reply by ●April 28, 20082008-04-28
On Apr 28, 5:08 pm, dbd <d...@ieee.org> wrote:> On Apr 28, 12:33 pm, Greg Heath <he...@alumni.brown.edu> wrote: > > Dft for Nonuniform Spacing > > > 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 > > 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?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
Reply by ●April 28, 20082008-04-28
On Apr 28, 7:23�pm, DigitalSignal <digitalsignal...@yahoo.com> wrote:> 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
Reply by ●April 29, 20082008-04-29
On Apr 28, 5:05 pm, Greg Heath <he...@alumni.brown.edu> wrote:> .... > 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. > ...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
Reply by ●April 29, 20082008-04-29
Greg Heath wrote:> Dft for Nonuniform Spacing > > 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. > > GregCheck Britz & Antonia, Meas. Sci. Technol. 7 (1996) 1042 and references therein. -- Dieter Britz (britz<at>chem.au.dk)
Reply by ●April 29, 20082008-04-29
On Apr 29, 12:27 am, dbd <d...@ieee.org> wrote:> On Apr 28, 5:05 pm, Greg Heath <he...@alumni.brown.edu> wrote: > > > .... > > 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. > > ... > > The question was and still is: > > > How closely does an evaluation of the spectrum information reproduce the > > original data? > > Jerry > > --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 A3Agree.> Hoping for help, > > Dale B. DalrympleThanks, again, for your previous help. I had been recommending AJJs dft code to others who needed code for nonuniform sampling. Greg
Reply by ●April 29, 20082008-04-29
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: > > > > .... > > > 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. > > > ... > > > The question was and still is: > > > > How closely does an evaluation of the spectrum information reproduce the > > > original data? > > > Jerry > > > -- > > 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. > > GregWhen 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
Reply by ●May 1, 20082008-05-01
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: > > > > > .... > > > > 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. > > > > ... > > > > The question was and still is: > > > > > How closely does an evaluation of the spectrum information reproduce the > > > > original data? > > > > Jerry > > > > -- > > > 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 > > 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.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?
