Prev: connecting to an existing matlab engine
Next: SimEvents randomising generator block seed in an model executable
From: Greg Heath on 25 Jul 2010 21:47
On May 27, 3:53 am, Greg Heath <he...(a)alumni.brown.edu> wrote:
> On May 25, 5:35 pm, Geoff <glmcd...(a)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 =     [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 thenon-> uniformdiscretefouriertransformto 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 thendftat
> > N = length(t_new);
> > f_new = (0:(N-1))/N;
>
> > % Use the modified AJ Johnson's implementation of thendftto evaluate
> > thendftat 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 inversefouriertransform(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)

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

 | 
Pages: 1
Prev: connecting to an existing matlab engine
Next: SimEvents randomising generator block seed in an model executable