Adjoint / inverse nufft

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

Answer 1

[Vilnis Liepins]

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.

[Victor Churchill]

Thanks for the response, Vilnis. The change you prescribe for the frequencies per the Nyquist frequency gave me exactly what I was looking for in your Figure 2. I'm not as familiar with this EDFT, but I'll check it out to see if it applies to my more nuanced problem. Thanks again!

[Paul]

@Vilnis Liepins

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]

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.

[Paul]

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]

t = 1×10

0 0.0047 0.1216 0.2784 0.4245 0.5434 0.6707 0.8259 0.8448 2.0000

Would the sampling frequency be considered as the maximum of the recriprocal of the samping intervals?

Fs = max(1./diff(t))

Fs = 211.9158

And then N is?

N = ceil(t(end)*Fs)

N = 424

and then the nufft f vector is?

f = (-ceil((N-1)/2):floor((N-1)/2))/N;

[Vilnis Liepins]

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.

[Paul]

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

Fs = 211.9158

Fs = 211.9158

N = ceil(t(end)*Fs)

N = 424

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].'

ans = 
   0.1367 + 0.0000i   2.3922 + 0.0044i
   0.5751 + 0.0000i   2.4081 + 0.0044i
   0.8913 + 0.0000i   2.7680 + 0.0039i
   0.2092 + 0.0000i   3.1195 + 0.0024i
   0.1853 + 0.0000i   3.2772 + 0.0006i
   0.1084 + 0.0000i   3.2719 - 0.0011i
   0.2197 + 0.0000i   3.1340 - 0.0027i
   0.9786 + 0.0000i   2.8005 - 0.0040i
   0.8117 + 0.0000i   2.7492 - 0.0041i
   0.1719 + 0.0000i  -0.2633 + 0.0044i

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.

[Paul]

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

Ts = 0.0100

Define the indices of the samples that correspond to the values the time vector based on this sampling period

nvec = round(t/Ts)

nvec = 1×10

0 10 25 30 33 49 53 58 70 100

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].'

ans = 
   0.5434 + 0.0000i   0.5434 - 0.0000i
   0.2784 + 0.0000i   0.2784 - 0.0000i
   0.4245 + 0.0000i   0.4245 - 0.0000i
   0.8448 + 0.0000i   0.8448 - 0.0000i
   0.0047 + 0.0000i   0.0047 - 0.0000i
   0.1216 + 0.0000i   0.1216 - 0.0000i
   0.6707 + 0.0000i   0.6707 - 0.0000i
   0.8259 + 0.0000i   0.8259 + 0.0000i
   0.1367 + 0.0000i   0.1367 + 0.0000i
   0.5751 + 0.0000i   0.5751 - 0.0000i

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)

dt = 1×9

0.0047 0.1169 0.1568 0.1461 0.1189 0.1273 0.1551 0.0189 1.1552

Ts = double(gcd(sym(dt)))

Ts = 1.1102e-16

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)

Ts = 0.0047

How close is t/Ts to integer valued?

format short e
[(t/Ts).' , (t/Ts).' - round(t/Ts).']

ans = 10×2

0 0 1.0000e+00 0 2.5762e+01 -2.3759e-01 5.8991e+01 -9.1400e-03 8.9962e+01 -3.8032e-02 1.1516e+02 1.5607e-01 1.4214e+02 1.4230e-01 1.7501e+02 1.1215e-02 1.7902e+02 2.1377e-02 4.2383e+02 -1.6848e-01

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].'

ans = 
   1.3671e-01 + 0.0000e+00i   1.5311e-01 + 6.0375e-03i
   5.7509e-01 + 0.0000e+00i   5.8186e-01 + 3.2108e-04i
   8.9132e-01 + 0.0000e+00i   8.9306e-01 + 5.0965e-03i
   2.0920e-01 + 0.0000e+00i   2.1546e-01 + 2.9315e-03i
   1.8533e-01 + 0.0000e+00i   1.9032e-01 + 7.1379e-04i
   1.0838e-01 + 0.0000e+00i   1.1169e-01 - 1.1398e-03i
   2.1970e-01 + 0.0000e+00i   2.1902e-01 - 3.5209e-03i
   9.7862e-01 + 0.0000e+00i   9.8172e-01 - 9.0553e-05i
   8.1168e-01 + 0.0000e+00i   8.1542e-01 - 2.8112e-04i
   1.7194e-01 + 0.0000e+00i   2.2665e-01 - 2.9570e-02i

figure
plot(t,x,t,real(x0),'-o')

[Vilnis Liepins]

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]

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.

Answer 2

[Matt J]

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

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

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

Answer 3

[Matt J]

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,

https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

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

[Victor Churchill]

Thanks for both of your answers, Matt J. While they indeed have expanded my understanding of the function, my priority here was inverting the nufft using the nufft function itself. Thanks again!