I am no expert on the non-uniform discrete fourier transform, but here
is the general goal of what I am working on. I have the input signal:
y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
missing samples at time 13 and 14
t =      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';

I want to reconstruct the two missing samples at time 13 and 14:
t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';

So the procedure I am trying, but not succeeding at is using the non-
uniform discrete fourier transform to evaluate the uniformly spaced
frequency domain, followed by a ifft() to reconstruct my new signal.
Here is the simple code I am trying:

---------CODE START:---------
% Signals
y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
t =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';

% Generate the frequencies I need to evaluate the ndft at
N = length(t_new);
f_new = (0:(N-1))/N;

% Use the modified AJ Johnson's implementation of the ndft to evaluate
the ndft at the needed frequency samples
y_tmp = y;
y_tmp(1) = 2*y_tmp(1);  % Scale, double first and last sample
y_tmp(end) = 2*y_tmp(end);
y_tmp = y_tmp/mean(diff(t));
[XFT,XLS,NMSEFT,NMSELS] = DFTgh1(y_tmp,t,f_new);

% Calculate in the inverse fourier transform (should be uniform
spacing in frequency domain and no ifftshift required)
y_new = ifft(XFT);
y_new_lsq = ifft(XLS);
------- CODE END---------

For some reason I can't get it to properly reconstruct the missing
samples. Instead it is placing those two missing samples somewhere
near a value of 0. I something wrong with my logic or understanding? I
know the code above does not work for even sample numbers, I just want
to get it working for the even case before considering the odd case.
(I would love to hear from the Greg Heath master on this topic as
well :) )


Appendix:
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

On May 26, 9:35=A0am, Geoff <glmcd...@gmail.com> wrote:
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =3D =A0 =A0 =A0[1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
>
> ---------CODE START:---------
> % Signals
> y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> % Generate the frequencies I need to evaluate the ndft at
> N =3D length(t_new);
> f_new =3D (0:(N-1))/N;
>
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp =3D y;
> y_tmp(1) =3D 2*y_tmp(1); =A0% Scale, double first and last sample
> y_tmp(end) =3D 2*y_tmp(end);
> y_tmp =3D y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new =3D ifft(XFT);
> y_new_lsq =3D ifft(XLS);
> ------- CODE END---------
>
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
>
> Appendix:
> function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)
>
> % function [XFT,XLS,NMSEFT,NMSELS] =3D 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) =3D sum(k=3D1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
> % =3D 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) =3D (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) =3D 0, then the MSE is unnormalized.
> %
> % NMSELS is the normalized mean-square-error of reconstucting
> % x from X using Least Squares. If mean(x) =3D 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 =3D
> % mean(diff(t))
>
> x =3D x(:); % Format 'x' into a column vector
> t =3D t(:); % Format 't' into a column vector
> f =3D f(:); % Format 'f' into a column vector
>
> N =3D length(x);
> if length(t)~=3D N
> =A0 =A0error('x and t do not have the same length')
> end;
>
> dt =3D diff(t); % asymmetric "dt"
> dts =3D 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
> meanx =3D sum(x.*dts)/sum(dts);
>
> df =3D diff(f); % asymmetric "df"
> dfs =3D 0.5*([df; 0]+[0; df]); % symmetric "df"
>
> W =3D exp(-2*pi*j * f*t');
>
> XFT =3D W * (x.*dts);
> XLS =3D (W'\x)./dfs;
>
> xFT =3D real(W'*(XFT.*dfs));
> xLS =3D real(W'*(XLS.*dfs));
>
> MSE0 =3D mse(abs(x-meanx));
> if MSE0 =3D=3D 0, MSE0 =3D 1, end;
>
> NMSEFT =3D mse(abs(x-xFT))/MSE0;
> NMSELS =3D mse(abs(x-xLS))/MSE0;
>
> return

What a palaver!!
Why not just use interp1?
And if its not smooth enough, use spline or pchip.

On May 25, 3:47=A0pm, TideMan <mul...@gmail.com> wrote:
> What a palaver!!
> Why not just use interp1?
> And if its not smooth enough, use spline or pchip.

Thanks, I will have a look at those functions more in-detail and see
if they will work for my application. My application is a bit special,
so it may not work as I need it to. I will post back with an update on
whether it works.

On May 26, 10:16=A0am, Geoff <glmcd...@gmail.com> wrote:
> On May 25, 3:47=A0pm, TideMan <mul...@gmail.com> wrote:
>
> > What a palaver!!
> > Why not just use interp1?
> > And if its not smooth enough, use spline or pchip.
>
> Thanks, I will have a look at those functions more in-detail and see
> if they will work for my application. My application is a bit special,
> so it may not work as I need it to. I will post back with an update on
> whether it works.

Yeah, well, I deal with data from tide gauges that have gaps.
These must be filled in order to do spectral or wavelet calculations.
One way to do it is to:
1. Put NaNs in the gaps.
2. Forecast the tide and subtract from the record to form a residual
3. Interpolate across the gaps using Rich Pawlowicz's elegant function
fixgaps (see below)
4. Add the forecast tide back in to the gap-filled residual.

function y=3Dfixgaps(x);
% FIXGAPS Linearly interpolates gaps in a time series
%      YOUT=3DFIXGAPS(YIN) linearly interpolates over NaN
%      in the input time series (may be complex), but ignores
%      trailing and leading NaN.
%

% R. Pawlowicz 6/Nov/99

y=3Dx;

bd=3Disnan(x);
gd=3Dfind(~bd);

bd([1:(min(gd)-1) (max(gd)+1):end])=3D0;


y(bd)=3Dinterp1(gd,x(gd),find(bd));

On May 25, 5:35=A0pm, Geoff <glmcd...@gmail.com> wrote:
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =3D =A0 =A0 =A0[1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
>
> ---------CODE START:---------
> % Signals
> y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> % Generate the frequencies I need to evaluate the ndft at
> N =3D length(t_new);
> f_new =3D (0:(N-1))/N;
>
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp =3D y;
> y_tmp(1) =3D 2*y_tmp(1); =A0% Scale, double first and last sample
> y_tmp(end) =3D 2*y_tmp(end);
> y_tmp =3D y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new =3D ifft(XFT);
> y_new_lsq =3D ifft(XLS);
> ------- CODE END---------
>
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
>
> Appendix:
> function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)
>
> % function [XFT,XLS,NMSEFT,NMSELS] =3D 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) =3D sum(k=3D1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
> % =3D 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) =3D (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) =3D 0, then the MSE is unnormalized.

CORRECTION:

% If mean(x) =3D x (i.e., x=3Dconstant) then the MSE is unnormalized.

> % NMSELS is the normalized mean-square-error of reconstucting
> % x from X using Least Squares. If mean(x) =3D 0, then the MSE
> % is unnormalized.

CORRECTION:

% If mean(x) =3D x (i.e., x=3Dconstant) 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 =3D
> % mean(diff(t))
>
> x =3D x(:); % Format 'x' into a column vector
> t =3D t(:); % Format 't' into a column vector
> f =3D f(:); % Format 'f' into a column vector
>
> N =3D length(x);
> if length(t)~=3D N
> =A0 =A0error('x and t do not have the same length')
> end;
>
> dt =3D diff(t); % asymmetric "dt"
> dts =3D 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
> meanx =3D sum(x.*dts)/sum(dts);
>
> df =3D diff(f); % asymmetric "df"
> dfs =3D 0.5*([df; 0]+[0; df]); % symmetric "df"
>
> W =3D exp(-2*pi*j * f*t');
>
> XFT =3D W * (x.*dts);
> XLS =3D (W'\x)./dfs;
>
> xFT =3D real(W'*(XFT.*dfs));
> xLS =3D real(W'*(XLS.*dfs));

CORRECTION:

xFT =3D real( W'*(XLS.*dfs) );
xLS =3D real( (W\XFT)./dts) );


> MSE0 =3D mse(abs(x-meanx));
> if MSE0 =3D=3D 0, MSE0 =3D 1, end;
>
> NMSEFT =3D mse(abs(x-xFT))/MSE0;
> NMSELS =3D mse(abs(x-xLS))/MSE0;
>
> return

The ifft is not appropriate for inverting the ndft
when dt is nonuniform. As indicated by the above
corrections, the correct inverse to XFT is xLS
obtained using least squares. In order to use
it for interpolation you would have to construct
a new W that depends on fnew and tnew.

Similarly for XLS and xFT. However, in this
case, perhaps ifft can be used.

This is the first time that I have thought of this
and I'm not positive it will work. I cannot get
back to thinking about this until next week.
Meanwhile, if you have the time, you might
try to pursue using XLS in ifft.

Greg

P.S. I'm tired and sleepy and cannot think
of when this would be more cost effective.
than other interpolation algorithms.

Geoff <glmcdona@gmail.com> wrote in message <a725007e-b73e-4ef7-b1ee-f3585c45bd41@p5g2000pri.googlegroups.com>...
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> 
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
> 
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
> 
> ---------CODE START:---------
> % Signals
> y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
> 
> % Generate the frequencies I need to evaluate the ndft at
> N = length(t_new);
> f_new = (0:(N-1))/N;
> 
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp = y;
> y_tmp(1) = 2*y_tmp(1);  % Scale, double first and last sample
> y_tmp(end) = 2*y_tmp(end);
> y_tmp = y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] = DFTgh1(y_tmp,t,f_new);
> 
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new = ifft(XFT);
> y_new_lsq = ifft(XLS);
> ------- CODE END---------
> 
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
> 
> 
> Appendix:
> 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
> 
Try 'edft' code available on fileexchange http://www.mathworks.com/matlabcentral/fileexchange/11020-extended-dft. 
Run command line: 
real(ifft(edft([1 2 3 4 5 6 7 8 9 10 11 12 NaN NaN 15 16 17]))) . 
It should return: 
ans =
  Columns 1 through 7 
    1.0000    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000
  Columns 8 through 14 
    8.0000    9.0000   10.0000   11.0000   12.0000   12.1150   13.0504
  Columns 15 through 17 
   15.0000   16.0000   17.0000
I suggest to increase DFT length twice:
real(ifft(edft([1 2 3 4 5 6 7 8 9 10 11 12 NaN NaN 15 16 17],2*17)))
ans =
  Columns 1 through 7 
    1.0000    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000
  Columns 8 through 14 
    8.0000    9.0000   10.0000   11.0000   12.0000   13.0010   14.0015
  Columns 15 through 21 
   15.0000   16.0000   17.0000   17.9634   18.7881   19.3011   19.2860
  Columns 22 through 28 
   18.5416   16.9515   14.5398   11.4925    8.1331    4.8557    2.0352
  Columns 29 through 34 
   -0.0565   -1.3002   -1.7353   -1.5202   -0.8680    0.0199
The first 17 samples is the reconstructed signal and samples 18..34 - the extapolation.

Thanks for all the replies. I will be taking a closer look at this later this week and will post back here.

Hi !

I am working on the same problem ! I just posted a simple algorithm which should work bur not working for me ! please have a look

Thanks

"Geoff McDonald" <glmcdona@gmail.com> wrote in message <hu10gj$gng$1@fred.mathworks.com>...
> Thanks for all the replies. I will be taking a closer look at this later this week and will post back here.

Also what i learnt form the guidance of people here ! LIke greag ! Walter. ! Tideman! and many other (thanks to all of them) you cant simply reconstruct the data you need to follow some regularization technique or need to do some thing Fourier domain ! to get uniform signal 

Geoff <glmcdona@gmail.com> wrote in message <a725007e-b73e-4ef7-b1ee-f3585c45bd41@p5g2000pri.googlegroups.com>...
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> 
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
> 
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
> 
> ---------CODE START:---------
> % Signals
> y =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
> 
> % Generate the frequencies I need to evaluate the ndft at
> N = length(t_new);
> f_new = (0:(N-1))/N;
> 
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp = y;
> y_tmp(1) = 2*y_tmp(1);  % Scale, double first and last sample
> y_tmp(end) = 2*y_tmp(end);
> y_tmp = y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] = DFTgh1(y_tmp,t,f_new);
> 
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new = ifft(XFT);
> y_new_lsq = ifft(XLS);
> ------- CODE END---------
> 
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
> 
> 
> Appendix:
> 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
> 

On May 27, 3:53=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On May 25, 5:35=A0pm, Geoff <glmcd...@gmail.com> wrote:
> > I am no expert on thenon-uniformdiscretefouriertransform, but here
> > is the general goal of what I am working on. I have the input signal:
> > y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> > missing samples at time 13 and 14
> > t =3D =A0 =A0 =A0[1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>
> > I want to reconstruct the two missing samples at time 13 and 14:
> > t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> > So the procedure I am trying, but not succeeding at is using thenon-> u=
niformdiscretefouriertransformto evaluate the uniformly spaced
> > frequency domain, followed by a ifft() to reconstruct my new signal.
> > Here is the simple code I am trying:
>
> > ---------CODE START:---------
> > % Signals
> > y =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> > t =3D =A0 =A0 [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> > t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>
> > % Generate the frequencies I need to evaluate thendftat
> > N =3D length(t_new);
> > f_new =3D (0:(N-1))/N;
>
> > % Use the modified AJ Johnson's implementation of thendftto evaluate
> > thendftat the needed frequency samples
> > y_tmp =3D y;
> > y_tmp(1) =3D 2*y_tmp(1); =A0% Scale, double first and last sample
> > y_tmp(end) =3D 2*y_tmp(end);
> > y_tmp =3D y_tmp/mean(diff(t));
> > [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>
> > % Calculate in the inversefouriertransform(should be uniform
> > spacing in frequency domain and no ifftshift required)
> > y_new =3D ifft(XFT);
> > y_new_lsq =3D ifft(XLS);
> > ------- CODE END---------
>
> > For some reason I can't get it to properly reconstruct the missing
> > samples. Instead it is placing those two missing samples somewhere
> > near a value of 0. I something wrong with my logic or understanding? I
> > know the code above does not work for even sample numbers, I just want
> > to get it working for the even case before considering the odd case.
> > (I would love to hear from the Greg Heath master on this topic as
> > well :) )
>
> > Appendix:
> > function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)

-----SNIP

>
> The ifft is not appropriate for inverting the ndft
> when dt is nonuniform. As indicated by the above
> corrections, the correct inverse to XFT is xLS
> obtained using least squares. In order to use
> it for interpolation you would have to construct
> a new W that depends on fnew and tnew.
>
> Similarly for XLS and xFT. However, in this
> case, perhaps ifft can be used.
>
> This is the first time that I have thought of this
> and I'm not positive it will work. I cannot get
> back to thinking about this until next week.
> Meanwhile, if you have the time, you might
> try to pursue using XLS in ifft.

-----SNIP

1. The corrected dftgh1 code is valid for
   periodic inputs y(end) =3D y(1). However,
   this leads to an illconditioned W with
   the first and last rows approximately
   equal.
   Replace dftgh1 with dftgh6 (below).
2. min(t) ~=3D 0. Therefore the corresponding
   input to dftgh6 should be t-min(t).
3. Least square spectra Xqr and Xpi (but not
   the DFT spectrum Xdft) are to be inverted
   using ifft.
4. The discontinuity abs(y(end)-y(1)) is so
   large that windowing is required before
   transforming and/or inverse transforming.
5. To interpolate in time, manually zeropad
   the spectra at the high frequency ends
   without losing conjugate symmetry.
6. Do not use the second argument in ifft(X,M)
   because it does not zeropad spectra correctly.
7. Techniques in other posts based on linearly
   interpolating over missing data before
   transforming should obviously work well for
   this linear function. Nevertheless, windowing
   is still required for good interpolation over
   the complete time interval.

Hope this helps.

Greg

function [Xdft,Xqr,Xpi,NMSEx] =3D dftgh6(x,t,f)

% function [Xdft,Xqr,Xpi,NMSEx] =3D dftgh6(x,t,f)
%
% Extension of AJ Johnson's dft to NONUNIFORM sampling
%
% Compute Xdft (Discrete Fourier Transform) at M frequencies
% given in f, given N samples x taken at times t:
%
% Xdft(f) =3D sum(k=3D1,N){dts(k)*x(k)*exp(-j*2*pi*f*t(k)*f)}
%         =3D W *(x.*dts)
% where dts is a symmetrized modification of diff(t).
%
% Compute least squares spectra Xqr and Xpi using
% using QR backslash and pseudoinverse least square
% inversion, respectively.
%
% NMSEx is the normalized mean-square-error of
% reconstucting x from the Xdft, Xqr and Xpi spectra.
% If meanx =3D x, i.e., x is constant, then the MSE is
% unnormalized.
%
% NMSEx(1) error from  Xqr  =3D>  xidftqr
% NMSEx(2) error from  Xpi  =3D>  xidftpi
% NMSEx(3) error from  Xdft =3D>  xiqrdft
% NMSEx(4) error from  Xdft =3D>  xipidft
%
% For comparison with MATLAB's Xfft when the spacing is
% uniform, multiply Xfft by dt0 =3D mean(diff(t))

x =3D x(:); % Format 'x' into a column vector
t =3D t(:); % Format 't' into a column vector
f =3D f(:); % Format 'f' into a column vector
N =3D length(x);
if length(t)~=3D N
   error('x and t do not have the same length')
end;

dt    =3D diff(t);                % asymmetric "dt"
% symmetric "dt"
dts   =3D [dt(1);(dt(1:end-1)+dt(2:end))/2;dt(end)];
meanx =3D sum(x.*dts)/sum(dts);

W    =3D exp(-2*pi*j * f*t');
Xdft =3D W * (x.*dts);

df   =3D diff(f);                 % asymmetric "df"
% symmetric  "df"
dfs  =3D [df(1);(df(1:end-1)+df(2:end))/2;df(end)];

% x =3D (1/N) *W' *(X.*dfs)       IDFT

Xqr  =3D ( W'\x )./dfs;
Xpi  =3D (pinv(W')*x)./dfs;

xidftqr =3D W' *(Xqr.*dfs);         % Fourier Reconstruction
xidftpi =3D W' *(Xpi.*dfs);         % Fourier Reconstruction

xiqrdft =3D ( W\Xdft )./dts ;       % QR LSE Reconstruction
xipidft =3D ( pinv(W)*Xdft )./dts ; % PI LSE Reconstruction

MSE0  =3D mse(abs(x-meanx));
if MSE0 =3D=3D 0, MSE0 =3D 1, end;

NMSEx(1)=3D mse(abs(x-xidftqr))/MSE0;
NMSEx(2)=3D mse(abs(x-xidftpi))/MSE0;
NMSEx(3)=3D mse(abs(x-xiqrdft))/MSE0;
NMSEx(4)=3D mse(abs(x-xipidft))/MSE0;

return

El martes 25 de mayo de 2010 18:35:28 UTC-3, Geoff  escribi=F3:
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =3D      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>=20
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
>=20
> ---------CODE START:---------
> % Signals
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> % Generate the frequencies I need to evaluate the ndft at
> N =3D length(t_new);
> f_new =3D (0:(N-1))/N;
>=20
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp =3D y;
> y_tmp(1) =3D 2*y_tmp(1);  % Scale, double first and last sample
> y_tmp(end) =3D 2*y_tmp(end);
> y_tmp =3D y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>=20
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new =3D ifft(XFT);
> y_new_lsq =3D ifft(XLS);
> ------- CODE END---------
>=20
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
>=20
>=20
> Appendix:
> function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)
>=20
> % function [XFT,XLS,NMSEFT,NMSELS] =3D 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) =3D sum(k=3D1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
> % =3D 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) =3D (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) =3D 0, then the MSE is unnormalized.
> %
> % NMSELS is the normalized mean-square-error of reconstucting
> % x from X using Least Squares. If mean(x) =3D 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 =3D
> % mean(diff(t))
>=20
> x =3D x(:); % Format 'x' into a column vector
> t =3D t(:); % Format 't' into a column vector
> f =3D f(:); % Format 'f' into a column vector
>=20
> N =3D length(x);
> if length(t)~=3D N
>    error('x and t do not have the same length')
> end;
>=20
> dt =3D diff(t); % asymmetric "dt"
> dts =3D 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
> meanx =3D sum(x.*dts)/sum(dts);
>=20
> df =3D diff(f); % asymmetric "df"
> dfs =3D 0.5*([df; 0]+[0; df]); % symmetric "df"
>=20
> W =3D exp(-2*pi*j * f*t');
>=20
> XFT =3D W * (x.*dts);
> XLS =3D (W'\x)./dfs;
>=20
> xFT =3D real(W'*(XFT.*dfs));
> xLS =3D real(W'*(XLS.*dfs));
>=20
> MSE0 =3D mse(abs(x-meanx));
> if MSE0 =3D=3D 0, MSE0 =3D 1, end;
>=20
> NMSEFT =3D mse(abs(x-xFT))/MSE0;
> NMSELS =3D mse(abs(x-xLS))/MSE0;
>=20
> return

testing

Hi Greg,

I have a question regarding with points 4 and 5 that you wrote in the last =
post.
What type of windowing is necessary? My idea is to do a Fourier interpolati=
on, take ndft and then inverse it to get original interpolated.
I am trying to zeropad manually, but i dont get quite good the idea.. is so=
mthing like this:

-code
L=3Dt(end);   % L=3D1800, last value for padding zeroes
t=3Dt-min(t); % now t(0)=3D0
% compute edft
[Xdft,Xqr,Xpi,NMSEx] =3D dftgh6(x,t,f);
%zero padding
% Xqrz =3D [Xqr(1:N/2); zeros(L-N,1); Xqr(N/2+1:N)];
% Xpiz =3D [Xpi(1:N/2); zeros(L-N,1); Xpi(N/2+1:N)];

Thanks in advance!

Alex

El martes 25 de mayo de 2010 18:35:28 UTC-3, Geoff  escribi=F3:
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =3D      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>=20
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
>=20
> ---------CODE START:---------
> % Signals
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> % Generate the frequencies I need to evaluate the ndft at
> N =3D length(t_new);
> f_new =3D (0:(N-1))/N;
>=20
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp =3D y;
> y_tmp(1) =3D 2*y_tmp(1);  % Scale, double first and last sample
> y_tmp(end) =3D 2*y_tmp(end);
> y_tmp =3D y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>=20
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new =3D ifft(XFT);
> y_new_lsq =3D ifft(XLS);
> ------- CODE END---------
>=20
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
>=20
>=20
> Appendix:
> function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)
>=20
> % function [XFT,XLS,NMSEFT,NMSELS] =3D 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) =3D sum(k=3D1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
> % =3D 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) =3D (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) =3D 0, then the MSE is unnormalized.
> %
> % NMSELS is the normalized mean-square-error of reconstucting
> % x from X using Least Squares. If mean(x) =3D 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 =3D
> % mean(diff(t))
>=20
> x =3D x(:); % Format 'x' into a column vector
> t =3D t(:); % Format 't' into a column vector
> f =3D f(:); % Format 'f' into a column vector
>=20
> N =3D length(x);
> if length(t)~=3D N
>    error('x and t do not have the same length')
> end;
>=20
> dt =3D diff(t); % asymmetric "dt"
> dts =3D 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
> meanx =3D sum(x.*dts)/sum(dts);
>=20
> df =3D diff(f); % asymmetric "df"
> dfs =3D 0.5*([df; 0]+[0; df]); % symmetric "df"
>=20
> W =3D exp(-2*pi*j * f*t');
>=20
> XFT =3D W * (x.*dts);
> XLS =3D (W'\x)./dfs;
>=20
> xFT =3D real(W'*(XFT.*dfs));
> xLS =3D real(W'*(XLS.*dfs));
>=20
> MSE0 =3D mse(abs(x-meanx));
> if MSE0 =3D=3D 0, MSE0 =3D 1, end;
>=20
> NMSEFT =3D mse(abs(x-xFT))/MSE0;
> NMSELS =3D mse(abs(x-xLS))/MSE0;
>=20
> return



El martes 25 de mayo de 2010 18:35:28 UTC-3, Geoff  escribi=F3:
> I am no expert on the non-uniform discrete fourier transform, but here
> is the general goal of what I am working on. I have the input signal:
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]'; % Straight line
> missing samples at time 13 and 14
> t =3D      [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
>=20
> I want to reconstruct the two missing samples at time 13 and 14:
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> So the procedure I am trying, but not succeeding at is using the non-
> uniform discrete fourier transform to evaluate the uniformly spaced
> frequency domain, followed by a ifft() to reconstruct my new signal.
> Here is the simple code I am trying:
>=20
> ---------CODE START:---------
> % Signals
> y =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t =3D     [1 2 3 4 5 6 7 8 9 10 11 12 15 16 17]';
> t_new =3D [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17]';
>=20
> % Generate the frequencies I need to evaluate the ndft at
> N =3D length(t_new);
> f_new =3D (0:(N-1))/N;
>=20
> % Use the modified AJ Johnson's implementation of the ndft to evaluate
> the ndft at the needed frequency samples
> y_tmp =3D y;
> y_tmp(1) =3D 2*y_tmp(1);  % Scale, double first and last sample
> y_tmp(end) =3D 2*y_tmp(end);
> y_tmp =3D y_tmp/mean(diff(t));
> [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(y_tmp,t,f_new);
>=20
> % Calculate in the inverse fourier transform (should be uniform
> spacing in frequency domain and no ifftshift required)
> y_new =3D ifft(XFT);
> y_new_lsq =3D ifft(XLS);
> ------- CODE END---------
>=20
> For some reason I can't get it to properly reconstruct the missing
> samples. Instead it is placing those two missing samples somewhere
> near a value of 0. I something wrong with my logic or understanding? I
> know the code above does not work for even sample numbers, I just want
> to get it working for the even case before considering the odd case.
> (I would love to hear from the Greg Heath master on this topic as
> well :) )
>=20
>=20
> Appendix:
> function [XFT,XLS,NMSEFT,NMSELS] =3D DFTgh1(x,t,f)
>=20
> % function [XFT,XLS,NMSEFT,NMSELS] =3D 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) =3D sum(k=3D1,N){ dts(k) *x(k) * exp(-2*pi*j*t(k)*f) }
> % =3D 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) =3D (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) =3D 0, then the MSE is unnormalized.
> %
> % NMSELS is the normalized mean-square-error of reconstucting
> % x from X using Least Squares. If mean(x) =3D 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 =3D
> % mean(diff(t))
>=20
> x =3D x(:); % Format 'x' into a column vector
> t =3D t(:); % Format 't' into a column vector
> f =3D f(:); % Format 'f' into a column vector
>=20
> N =3D length(x);
> if length(t)~=3D N
>    error('x and t do not have the same length')
> end;
>=20
> dt =3D diff(t); % asymmetric "dt"
> dts =3D 0.5*([dt; 0]+[0; dt]); % symmetric "dt"
> meanx =3D sum(x.*dts)/sum(dts);
>=20
> df =3D diff(f); % asymmetric "df"
> dfs =3D 0.5*([df; 0]+[0; df]); % symmetric "df"
>=20
> W =3D exp(-2*pi*j * f*t');
>=20
> XFT =3D W * (x.*dts);
> XLS =3D (W'\x)./dfs;
>=20
> xFT =3D real(W'*(XFT.*dfs));
> xLS =3D real(W'*(XLS.*dfs));
>=20
> MSE0 =3D mse(abs(x-meanx));
> if MSE0 =3D=3D 0, MSE0 =3D 1, end;
>=20
> NMSEFT =3D mse(abs(x-xFT))/MSE0;
> NMSELS =3D mse(abs(x-xLS))/MSE0;
>=20
> return

[testing the new google.groups iface]

Hi Greg,

I have a question regarding with points 4 and 5 that you wrote in the last =
post.
What type of windowing is necessary? My idea is to do a Fourier interpolati=
on, take ndft and then inverse it to get original interpolated.
I am trying to zeropad manually, but i dont get quite good the idea.. is so=
mthing like this:

----code.m
....
L=3Dt(end);   % L=3D1800, last value for padding zeroes
t=3Dt-min(t); % now t(0)=3D0
% compute edft
[Xdft,Xqr,Xpi,NMSEx] =3D dftgh6(x,t,f);
%zero padding
% Xqrz =3D [Xqr(1:N/2); zeros(L-N,1); Xqr(N/2+1:N)];
% Xpiz =3D [Xpi(1:N/2); zeros(L-N,1); Xpi(N/2+1:N)];
....

Thanks in advance!

Sorry if I resend this message, I tried a couple times unsuccessfully...

Alex

On Feb 29, 10:12=A0pm, Alex Santa-Cruz <websan...@gmail.com> wrote:

---SNIP

Thanks for emailing me your code and data.

I rewrote your program to use dftgh6 properly. Unfortunately, the
bottom line is
that iffting the zeropadded Xpi (Xqr doesn't work at all) dft spectrum
just produces
an interpolated version of the original signal with zeros
(approximately) filling the
gaps.

This makes sense since it is the same dft/fft spectrum one would get
if the original
signal gaps were to be filled with zeros.

So, this technique is OK for producing an interpolated uniformly
spaced version
of a nonuniformly spaced signal, provided the spaces between the
original data
points are not too large.

My memory is not good these days however, I think that I was able to
use this
to fill the gap in the 15 point straight line [ 1:12, 15:17 ] of the
OP. I "remember"
shifting the time to start at 0 and reducing the Gibbs phenomenom by
replacing
the 17 with (17-1)/2. Unfortunately, those results are on a previous
computer that
died.

Other posters mentioned
1. Using linear interpolation to fill gaps (obviously works for the
OPs line) before
   using the dft
2. Using edft from the MATLAB File-Exchange.

The latter uses some sort of spectral prediction method for
interpolation and
extrapolation.

I think this may be worth investigating. Circa 1978 a colleague of
mine constructed
a gap filler based on spectral linear prediction. Not sure if there
were any unclassified
publications. Search on "Steven/Stephen Bowling".

My best advice is to down load edft and try to duplicate the straight
line example in
the thread as well as any examples in the edft package.

Hope this helps.

Greg