Prev: Specifying non-linear axis
Next: Plot random markers
From: Vilnis Liepins on 27 May 2010 10:18
Geoff <glmcdona(a)gmail.com> wrote in message <a725007e-b73e-4ef7-b1ee-f3585c45bd41(a)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.
From: Geoff McDonald on 31 May 2010 14:52
Thanks for all the replies. I will be taking a closer look at this later this week and will post back here.
From: kk KKsingh on 31 May 2010 14:59
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(a)gmail.com> wrote in message <hu10gj$gng$1(a)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.
From: kk KKsingh on 31 May 2010 15:52
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(a)gmail.com> wrote in message <a725007e-b73e-4ef7-b1ee-f3585c45bd41(a)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
>
 | 
Pages: 1
Prev: Specifying non-linear axis
Next: Plot random markers