Adjoint / inverse nufft
Show older comments
Here's a simple example of the behavior based on the documentation for nufft:
t = [0:300 500.5:700.5];
S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t);
X = abs(S + rand(size(t))); % true signal
n = length(t);
f = (0:n-1)/n;
Y = nufft(X,t,f);
X0 = real(nufft(Y,f,-t))/n; % reconstructed signal
figure,plot(t,X); hold on; plot(t,X0)
How can I compensate for the nonuniformity in t to get X0 to match X?
One answer to this post (link) mentions a "density compensation matrix", but no details and there are no other outputs from the nufft function. I assume this has something to do with the relationship (interpolation or whatever is going on behind the scenes) of the non-uniform to uniform grid.
This plot shows that the difference indeed differs over t.
figure,plot(t,abs(X-X0))
Accepted Answer
Vilnis Liepins
on 13 Oct 2023
1 vote
Hi Victor,
You can improve the inverse of NUFFT if take into account that the length of signal in time is 700.5 sec and choose, e.g,
n=701; You should also select frequencies that meet the Nyquist limit of 0.5 Hz, f=(-ceil((n-1)/2):floor((n-1)/2))/n;
However, the inverse of NUFFT will not match to X and if you apply a jitter to t = t + rand(1,length(t))/2; then the differences will get more visible.
If you want to get exactly the same signal back, I recommend trying the EDFT code in fileexchange https://se.mathworks.com/matlabcentral/fileexchange/11020-extended-dft , then you get picture
Moreover, you can apply Matlab IFFT to the result of EDFT and get the signal resampled on regular grid and the gap filled with interpolated values
My code
t = [0:300 500.5:700.5];
% t = t + rand(1,length(t))/2; % add jitter to t
S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t);
X = abs(S + rand(size(t)));
%n = length(t);
n=701; % The last sample taken at 700.5 sec
%f = (0:n-1)/n;
f=(-ceil((n-1)/2):floor((n-1)/2))/n; % Nyquist frequency = 0.5 Hz
% NUFFT
Y_NUFFT = nufft(X,t,f);
S_NUFFT=Y_NUFFT/length(X); % Power Spectrum
% EDFT
[Y_EDFT,S_EDFT,st]=edft(X,f,[],[],t);
figure(1) %Plot Power Spectrums
plot(f,abs(S_NUFFT)), hold on, plot(f,abs(S_EDFT)), hold off
X0 = real(nufft(Y_NUFFT,f,-t))/n; % reconstructed signal NUFFT
figure(2),plot(t,X), hold on, plot(t,X0), hold off
X1 = real(iedft(Y_EDFT,f,t)); % reconstructed signal IEDFT
figure(3),plot(t,X), hold on, plot(t,X1), hold off
X2 = real(ifft(ifftshift(Y_EDFT))); % reconstructed signal by IFFT on uniform grid
figure(4),plot(t,X), hold on, plot(0:length(f)-1,X2), hold off
I hope it helps. Don't hesitate to ask if you have any questions.
9 Comments
Victor Churchill
on 13 Oct 2023
Paul
on 14 Oct 2023
Can you explain why f has to be defined as
f=(-ceil((n-1)/2):floor((n-1)/2))/n;
instead of
f = (0:n-1)/n;
to get the desired x0
Vilnis Liepins
on 14 Oct 2023
Edited: Vilnis Liepins
on 14 Oct 2023
Short answer - because of signal is sampled nonuniformly and we would like to calculate the Fourier Transform defined for positive and negative frequencies.
If signal sampled in uniform way you can draw a sine wave with a frequency above Nyquist (> 0.5 Hz) and a negative one (< 0 Hz) througt exactly the same sample points. You can apply 'fftshift' to FFT output and the spectrum above 0.5 Hz will be shifted into negative frequency range - this works perfectly.
All of the above is not true if the sampling points are taken nonuniformly - NUFFT on frequencies above Nyquist are not the same as at negative frequencies. That's why using f=(-ceil((n-1)/2):floor((n-1)/2))/n; preferable.
Ah, I think I understand. I have trouble getting my mind wrapped around the NUDFT as a concept.
Follow up question. If the sampling points are not uniform, what do we use for the sample period for purpose of determining the sampling frequency and the Nyquist frequency?
For example, suppose I have a finite duration sequence derived by sampling a continuous time signal nonuniformly over 0 <= t < = 2
rng(100)
t = [0 sort(rand(1,8)) 2]
Would the sampling frequency be considered as the maximum of the recriprocal of the samping intervals?
Fs = max(1./diff(t))
And then N is?
N = ceil(t(end)*Fs)
and then the nufft f vector is?
f = (-ceil((N-1)/2):floor((N-1)/2))/N;
Vilnis Liepins
on 15 Oct 2023
Hi Paul,
Sampling frequency (Fs) and Nyquist frequency (Fs/2) could be determined from Mean Sampling Period (Ts).
This works the same as in uniform case, e.g., let resample your signal uniformly:
t = [0:0.1:0.8 2.0]
t =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 2.0000
What is the sampling period T? Obviously 0.1 sec and Fs = 1/T = 10 Hz. The gap in 1 sec and the last sample at t = 2 sec is excluded from the estimation.
Back to your example, in nonuniform case we calculate: mean(diff(t(1:9))) = 0.1056
and for simplicity we choose the value Ts=0.1 sec - the same as in uniform case. Then Fs = 1/Ts = 10 Hz and the frequency set could be generated as
f = Fs*(-ceil((N-1)/2):floor((N-1)/2))/N; and this covers one Nyquist zone.
I fact, nonuniform sampling allows you to estimate signal spectrum also beyond one Nyquist zone. For example, you can generate frequencies to cover 2 Nyquist zones:
Nq_zone=2;
f = Fs*(-ceil(Nq_zone*(N-1)/2):floor(Nq_zone*(N-1)/2))/N;
How to select the number of frequencies N? Let's continue similarly to the Matlab example. Since the last sample is taken after 2 seconds, a good choice would be a value of N equal to or greater than the number of samples taken in 0 <= t < = 2 with Ts = 0.1 seconds. So we get N >= 21. Your calculated value of N = 424 is fine too.
If N = 424, then let's continue the example with that
rng(100)
t = [0 sort(rand(1,8)) 2];
Fs = max(1./diff(t))
Fs = 211.9158
N = ceil(t(end)*Fs)
f = (-ceil((N-1)/2):floor((N-1)/2))/N;
Generate signal samples at the time points.
x = rand(1,numel(t));
Now, I'll take the nufft of x and invert it, same as shown above
x0 = nufft(nufft(x,t,f),f,-t)/N;
And compare
[x;x0].'
The reconstructed signal doesn't look anything like the original, even neglecting the small imaginary part, which isn't so small
figure
plot(t,x,'-o',t,real(x0),'-o')
What is the difference between this example and the @Victor Churchill's original problem where this inversion worked so well?
Also, any links to good references on NUDFT (not Wikipedia) would be appreciated if you have any.
I think the answer to my own question is the nufft will be invertible if the sampling period is chosen based on the GCD of the delta-T's in the time vector.
Define a simple time vector and data samples that go with it.
t = [0, 0.1, 0.25, 0.3, 0.33, 0.49, 0.53, 0.58, 0.7, 1];
rng(100)
x = rand(1,numel(t));
Define the sampling period as the GCD of the delta-T's.
dt = diff(t);
Ts = double(gcd(sym(dt)))
Define the indices of the samples that correspond to the values the time vector based on this sampling period
nvec = round(t/Ts)
The total number of samples that cover the time interval with a sampling period of Ts is
N = nvec(end) + 1;
Create a data vector that is zero filled between the samples in x
xexpanded = zeros(1,N);
xexpanded(nvec+1) = x;
Define the frequency vector in the usual way.
f = (0:N-1)/N/Ts;
Now take the NUDFT of the nonuniformly spaced samples and the DFT of the zero filled samples.These should be identical because the zeros in xexpanded make no contribution to the DFT of xexpanded, and the NUDFT of x, with this frequency vector, is the therefore the same as the DFT of xexpanded
Xnufft = nufft(x,t,f);
Xexpfft = fft(xexpanded);
Verfy NUDFT of x and DFT of xexpanded are the same.
figure
plot(f,abs(Xnufft),'-o',f,abs(Xexpfft))
figure
plot(f,angle(Xnufft),'-o',f,angle(Xexpfft))
Because of the invertability of the DFT, and that the DFT and NUDFT are thte same for this construction, it stands to reason that the NUDFT should be invertible as well
x0 = nufft(Xnufft,f,-t)/N;
[x;x0].'
figure
plot(t,x,t,real(x0),'o')
So, it looks like the NUDFT is exactly invertible if the frequency vector and N are chosen based on a sampling period that is the GCD of the individual sampling intervals.
In practice, that's not, in general, going to be workable. Returning to my original example.
rng(100)
t = [0 sort(rand(1,8)) 2];
x = rand(1,numel(t));
dt = diff(t)
Ts = double(gcd(sym(dt)))
Ts computed this way is less than eps, so not very useful. I'm thinking that the next best way to pick Ts would such that t/Ts are "close" to integer valued. Don't know how close is close enough. In this problem, dt(1) is two orders of magnitude smaller than the others, so just use that
Ts = dt(1)
How close is t/Ts to integer valued?
format short e
[(t/Ts).' , (t/Ts).' - round(t/Ts).']
Don't know if that's close enough, but let's assume it is.
Proceed as above
N = ceil(t(end)/Ts) + 1;
f = (0:N-1)/N/Ts;
x0 = nufft(nufft(x,t,f),f,-t)/N;
[x;x0].'
figure
plot(t,x,t,real(x0),'-o')
Vilnis Liepins
on 16 Oct 2023
Edited: Vilnis Liepins
on 16 Oct 2023
The answer on your comment https://se.mathworks.com/matlabcentral/answers/2032659-adjoint-inverse-nufft#comment_2923034
In my previous comment, i suggest to calculate Sampling Frequency based on Mean Sampling Period and get the value Fs=10 Hz for your data. The second point you miss is that frequencies should be calculated as f = Fs*(-ceil((N-1)/2):floor((N-1)/2))/N; and the picture i got looks much better
rng(100)
t = [0 sort(rand(1,8)) 2];
Fs=10;
N=424;
f = Fs*(-ceil((N-1)/2):floor((N-1)/2))/N;
x = rand(1,numel(t));
x0 = nufft(nufft(x,t,f),f,-t)/N;
figure
plot(t,x,'-o',t,real(x0),'-o')
Paul
on 16 Oct 2023
Ah yes. In this comment there was no explicit use of Fs in the expression for f, but now I see that Fs is used in this comment. In the former Fs=1 becasue of the nature of the data in the problem posted by @Victor Churchill. Thanks for pointing that out.
More Answers (2)
If these are your actual data sizes, an optimization approach seems to work not too badly:
t = [0:300 500.5:700.5];
S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t);
X = abs(S + rand(size(t))); % true signal
n = length(t);
f = (0:n-1)/n;
Y = nufft(X,t,f);
X0=fsolve(@(x)resFunc(x,t,f,Y), rand(size(X)) );
plot(t,X,'--b',t,X0,'xr')
function r=resFunc(x,t,f,Y)
r=Y-nufft(x,t,f);
r=[real(r);imag(r)];
end
If your t,f sampling is going to be reused, it may be gainful to use an algebraic inversion with the help of this FEX download,
Note that the C matrix depends only on t and f and can be re-used.
t = [0:300 500.5:700.5];
S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t);
X = abs(S + rand(size(t))); % true signal
n = length(t);
f = (0:n-1)/n;
Y = nufft(X,t,f);
C=func2mat(@(z) nufft(z,t,f), X);
C=[real(C);imag(C)];
d=[real(Y(:));imag(Y(:))];
X0=C\d;
plot(t,X,'-bx',t,X0,'-r')
1 Comment
Victor Churchill
on 13 Oct 2023
Categories
Find more on Transforms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
