ndft (non-uniform discrete fourier transform) and reconstructing missing samples
Geoff
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
TideMan
What a palaver!!
Why not just use interp1?
And if its not smooth enough, use spline or pchip.
Geoff
> 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.
TideMan
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=fixgaps(x);
% FIXGAPS Linearly interpolates gaps in a time series
% YOUT=FIXGAPS(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=x;
bd=isnan(x);
gd=find(~bd);
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
y(bd)=interp1(gd,x(gd),find(bd));
Greg Heath
CORRECTION:
% If mean(x) = x (i.e., x=constant) 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.
CORRECTION:
% If mean(x) = x (i.e., x=constant) 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));
CORRECTION:
xFT = real( W'*(XLS.*dfs) );
xLS = real( (W\XFT)./dts) );
> MSE0 = mse(abs(x-meanx));
> if MSE0 == 0, MSE0 = 1, end;
>
> NMSEFT = mse(abs(x-xFT))/MSE0;
> NMSELS = 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.
Vilnis Liepins
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.
Geoff McDonald
kk KKsingh
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" <glmc...@gmail.com> wrote in message <hu10gj$gng$1...@fred.mathworks.com>...
kk KKsingh
Geoff <glmc...@gmail.com> wrote in message <a725007e-b73e-4ef7...@p5g2000pri.googlegroups.com>...
Greg Heath
> On May 25, 5:35 pm, Geoff <glmcd...@gmail.com> wrote:
> > 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]';
>
> > 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 thendftat
> > N = length(t_new);
> > f_new = (0:(N-1))/N;
>
> > % Use the modified AJ Johnson's implementation of thendftto evaluate
> > 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);
>
> > 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)
-----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) = 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) ~= 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] = dftgh6(x,t,f)
% function [Xdft,Xqr,Xpi,NMSEx] = 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) = sum(k=1,N){dts(k)*x(k)*exp(-j*2*pi*f*t(k)*f)}
% = 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 = x, i.e., x is constant, then the MSE is
% unnormalized.
%
% NMSEx(1) error from Xqr => xidftqr
% NMSEx(2) error from Xpi => xidftpi
% NMSEx(3) error from Xdft => xiqrdft
% NMSEx(4) error from Xdft => xipidft
%
% For comparison with MATLAB's Xfft when the spacing is
% uniform, multiply Xfft 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"
% symmetric "dt"
dts = [dt(1);(dt(1:end-1)+dt(2:end))/2;dt(end)];
meanx = sum(x.*dts)/sum(dts);
W = exp(-2*pi*j * f*t');
Xdft = W * (x.*dts);
df = diff(f); % asymmetric "df"
% symmetric "df"
dfs = [df(1);(df(1:end-1)+df(2:end))/2;df(end)];
% x = (1/N) *W' *(X.*dfs) IDFT
Xqr = ( W'\x )./dfs;
Xpi = (pinv(W')*x)./dfs;
xidftqr = W' *(Xqr.*dfs); % Fourier Reconstruction
xidftpi = W' *(Xpi.*dfs); % Fourier Reconstruction
xiqrdft = ( W\Xdft )./dts ; % QR LSE Reconstruction
xipidft = ( pinv(W)*Xdft )./dts ; % PI LSE Reconstruction
MSE0 = mse(abs(x-meanx));
if MSE0 == 0, MSE0 = 1, end;
NMSEx(1)= mse(abs(x-xidftqr))/MSE0;
NMSEx(2)= mse(abs(x-xidftpi))/MSE0;
NMSEx(3)= mse(abs(x-xiqrdft))/MSE0;
NMSEx(4)= mse(abs(x-xipidft))/MSE0;
return
Alex Santa-Cruz
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 interpolation, 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 somthing like this:
-code
L=t(end); % L=1800, last value for padding zeroes
t=t-min(t); % now t(0)=0
% compute edft
[Xdft,Xqr,Xpi,NMSEx] = dftgh6(x,t,f);
%zero padding
% Xqrz = [Xqr(1:N/2); zeros(L-N,1); Xqr(N/2+1:N)];
% Xpiz = [Xpi(1:N/2); zeros(L-N,1); Xpi(N/2+1:N)];
Thanks in advance!
Alex
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 interpolation, 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 somthing like this:
----code.m
...
L=t(end); % L=1800, last value for padding zeroes
t=t-min(t); % now t(0)=0
% compute edft
[Xdft,Xqr,Xpi,NMSEx] = dftgh6(x,t,f);
%zero padding
% Xqrz = [Xqr(1:N/2); zeros(L-N,1); Xqr(N/2+1:N)];
% Xpiz = [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
Greg Heath
---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