Example python nfft fourier transform - Issues with signal reconstruction normalization

Question

I wrote a full working example for both nfft, and scipy.fft. In both cases I start with a simple 1D sinusoidal signal with a little noise, take the fourier transform, and then go backwards and reconstruct the original signal.

Here is my code as clean and readable as I could manage:

import numpy
import nfft
import scipy
import scipy.fft
import matplotlib.pyplot as plt

if True: #<--- Ensure non-global namespace
    #Define signal:
    x = -0.5 + numpy.random.rand(1000)
    #x = numpy.linspace(-.5, .5, 1000) #--> in case we want to run uniform time domain
    f = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( 1000 ) #Add some  'y' randomness to the sample

    #prepare wavenumbers for transform:
    N = len(x)
    k = - N // 2 + numpy.arange(N) 
    #print ('k', k) #---> Uniform Steps [-500, -499, ...0..., 499,500]

    f_k = nfft.nfft_adjoint(x, f, len(k), truncated=False )

    #plot transform
    plt.figure()
    plt.plot(k, f_k.real, label='real')
    plt.plot(k, f_k.imag, label='imag')
    plt.legend()

    #Reconstruct the original signal with nfft
    f_recon = nfft.nfft( x, f_k ) / 2000

    #Plot original vs reconstructed
    plt.figure()
    plt.title('nfft')
    plt.scatter(x, f, label='f(x)')
    plt.scatter(x, f_recon, label='f_recon(x)', marker='+')
    plt.legend()

if True: #<--- Ensure non-global namespace
    #Define signal:
    x = numpy.linspace(-.5, .5, 1000)
    f = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( 1000 ) #Add some 'y' randomness to the sample

    #prepare wavenumbers for transform:
    N = len(x)
    TimeSpacing = x[1] - x[0]
    k = scipy.fft.fftfreq(N, TimeSpacing)
    #print ('k', k) #---> Confusing steps: [0,1,...500,-500,-499,...-1]
    
    f_k = scipy.fft.fft(f)

    #plot transform
    plt.figure()
    plt.plot(k, f_k.real, label='real')
    plt.plot(k, f_k.imag, label='imag')
    plt.legend()

    #Reconstruct the original signal with scipy.fft
    f_recon = scipy.fft.ifft(f_k)

    #Plot original vs reconstructed
    plt.figure()
    plt.title('scipy.fft')
    plt.scatter(x, f, label='f(x)')
    plt.scatter(x, f_recon, label='f_recon(x)', marker='+')
    plt.legend()

plt.show()

Here are the relevant generated plots:

The nfft reconstruction seems to fail at normalizing. I arbitrarily divided the magnitudes by 2000 just to get them to plot well. What is the correct normalization constant?

The nfft also seems to not reproduce the original points. Even if I got hte normalization constant correct, there is no way I would get the original points back here.

What did I do wrong, and how do I fix it?

Answer 1

The above mentioned package does not implement a inverse nfft

The ndft is f_hat @ np.exp(-2j * np.pi * x * k[:, None]) The ndft_adjoint is f @ np.exp(2j * np.pi * k * x[:, None])

Let k = -N//2 + np.arange(N), and A = np.exp(-2j * np.pi * k * k[:, None])

A @ np.conj(A) = N * np.eye(N) (checked numerically)

Thus, for random x the adjoint transformation is equals to the inverse transform. The given reference paper provides a few options, I implemented Algorithm 1 CGNE, from page 9

import numpy as np # I have the habit to use np
def nfft_inverse(x, y, N, w = 1, L=100):
    f = np.zeros(N, dtype=np.complex128);
    r = y - nfft.nfft(x, f);
    p = nfft.nfft_adjoint(x, r, N);
    r_norm = np.sum(abs(r)**2 * w)
    for l in range(L):
        p_norm = np.sum(abs(p)**2 * w);
        alpha = r_norm / p_norm
        f += alpha * w * p;
        r = y - nfft.nfft(x, f)
        r_norm_2 = np.sum(abs(r)**2 * w)
        beta = r_norm_2 / r_norm
        p = beta * p + nfft.nfft_adjoint(x, w * r, N)
        r_norm = r_norm_2;
        #print(l, r_norm)
    return f;

The algorithm converges slowly and poorly

    plt.figure(figsize=(14, 7))
    plt.title('inverse nfft error histogram')
    #plt.scatter(x, f_hat, label='f(x)')
    h_hat = nfft_inverse(x, f, N, L = 1)
    plt.hist(f_hat - numpy.real(h_hat), bins=30, label='1 iteration')
    h_hat = nfft_inverse(x, f, N, L = 10)
    plt.hist(f_hat - numpy.real(h_hat), bins=30, label='10 iterations')
    h_hat = nfft_inverse(x, f, N, L = 1000)
    plt.hist(f_hat - numpy.real(h_hat), bins=30, label='1000 iterations')
    plt.xlabel('error')
    plt.ylabel('occurrencies')
    plt.legend()

I tried also to use scipy minimization, to minimize the residual ||nfft(x, f) - y||**2 explicitly

import numpy as np # the habit
import scipy.optimize
def nfft_gradient_descent(x, y, N, L=10, tol=1e-8, method='CG'):
    '''
    compute $min || A @ f - y ||**2 via gradient descent
    the gradient is
    
    `A^H @ (A @ f - y)`
    
    Multiply by A using nfft.nfft
    
    '''
    def cost(fpack):
        f = fpack[0::2] + 1j * fpack[1::2]
        u = np.sum(np.abs(nfft.nfft(x, f) - y)**2)
        return u
    def grad(fpack):
        f = fpack[0::2] + 1j * fpack[1::2]
        r = nfft.nfft(x, f) - y
        u = nfft.nfft_adjoint(x, r, N)
        return np.stack([np.real(u), np.imag(u)], axis=1).reshape(-1)
    
    x0 = np.zeros([N, 2])
    result = scipy.optimize.minimize(cost, x0=x0, jac=grad, tol=tol, method=method, options={'maxiter': L, 'disp': True})
    return result.x[0::2] + 1j * result.x[1::2];

The results look similar, you can try by different methods or parameters by your self if you want. But I believe the transformation is ill conditioned because transformed residual is considerably reduced but the residual on the reconstructed values is large.

Edit 1

Is it basically true that you found that the there is not a true inverse to the algorithm then? I cannot obtain my original points? x != nfft(nfft_adjoint(x))

Please check the section 2.3 of the reference paper

Numerical comparison

Cris Luengo answer metioned another possibility, that is, instead of reconstructing f at the points x, you may reconstruct a resampled version at equidistant points using the ifft. So you already have three options, and I will do a quick comparison. Bear in mind that the plot shown there is based on a NFFT computed in 16k samples, while here I am using 1k samples.

Since the FFT method uses different points we cannot compare to the original signal, what I will do is to compare with the harmonic function without the noise. The variance of the noise is 0.01, so an exact reconstruction would lead to this mean squared error.

N = 1024
x = -0.5 + numpy.random.rand(N)
f_hat = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( N ) #Add some  'y' randomness to the sample

k = - N // 2 + numpy.arange(N)
f = nfft.nfft(x, f_hat)

print('nfft_inverse')
h_hat = nfft_inverse(x, f, len(x), L = 10)
print('10 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_inverse(x, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_inverse(x, f, len(x), L = 1000)
print('1000 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
print('nfft_gradient_descent')
h_hat = nfft_gradient_descent(x, f, len(x), L = 10)
print('10 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 1000)
print('1000 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))


#Reconstruct the original at a spaced grid based on nfft result using ifft
f_recon = - numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k))) / (N / N2)
x_recon = k / N; 
print('using IFFT: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x_recon) - numpy.real(f_recon))**2))

Results:

nfft_inverse
10 iterations:  0.0798988590351581
100 iterations:  0.05136853850272318
1000 iterations:  0.037316315280700896
nfft_gradient_descent
10 iterations:  0.08832834348902704
100 iterations:  0.05901599049633016
1000 iterations:  0.043921864589484
using IFFT:  0.49044932854606377

Another way to view is

plt.plot(numpy.sin(10 * 2 * numpy.pi * x_recon), numpy.real(f_recon), '.', label='ifft')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x), numpy.real(nfft_gradient_descent(x, f, len(x), L = 5)), '.', label='gradient descent L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x), numpy.real(nfft_inverse(x, f, len(x), L = 5)), '.', label='nfft_inverse L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x), np.real(f_hat), '.', label='original')
plt.legend()

Even though the IFFT matrix is better conditioned, it gives results in a worse reconstruction of the signal. Also from the last plot it becomes more noticeable that there is a slight attenuation. Probably due to a systematic energy leak to the imaginary part (an error in my code??). Just a quick test, multiplying it by 1.3 gives better results

Answer 2

Bob has already posted an excellent answer, this is just to complement with some details I hope are instructive.

First, compare the two plots for the computed frequency components. Note how the one for the NFFT is much more noisy than the one for the regular FFT. You are estimating these frequency components for a sampled signal from noisy samples, in one case the samples are regularly spaced, in the other they are randomly spaced. It is a known result that regular sampling is more efficient than random sampling (efficient meaning that you need fewer samples to get the same amount of information). Thus, it is expected that the random sampling yield the more noisy result.

We can compute the "normal" inverse FFT from the frequency components estimated by the NFFT:

f_recon = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k)))
x_recon = numpy.linspace(-.5, .5, N)

I used ifftshift because NFFT defines the k to go from -N/2 to N/2-1, whereas the FFT defines it to go from 0 to N-1. ifftshift swaps the two halves of the signal to turn the first into the second (k from N/2 to N-1 is equal to -N/2 to -1). I also used fftshift on the result of the IFFT, because the same thing applies to the time axis, it shifts the origin from the first sample to the middle of the sequence.

Note how noisy f_recon is. This is explained by the poor estimates of f_k we could make with the non-uniformed sampled signal. There is also a sign error, we could already observe this error when we compared the two estimates of f_k. This comes from the adjoint NFFT having the same sign in the exponent as the inverse DFT, which actually means that f_recon is flipped w.r.t. x.

If we increase the number of random samples we can obtain a better estimate:

import numpy
import nfft
import matplotlib.pyplot as plt

#Define signal:
N = 1024 * 16  # power of two for speed
x = -0.5 + numpy.random.rand(N)
f = numpy.sin(10 * 2 * numpy.pi * x) + .1 * numpy.random.randn(N) # Add some  'y' randomness to the sample

#prepare wavenumbers for transform:
k = - N // 2 + numpy.arange(N) 

N2 = 1024
f_k = nfft.nfft_adjoint(x, f, N2, truncated=False)

#Reconstruct the original signal with nfft
#   (note the minus sign to flip the signal, in reality we should flip x)
f_recon = - numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k))) / (N / N2)
x_recon = numpy.linspace(-.5, .5, N2, endpoint=False)

#Plot original vs reconstructed
plt.figure()
plt.title('nfft')
plt.scatter(x[:N2], f[:N2], label='f(x)') # don't plot all samples, there's too many
plt.scatter(x_recon, f_recon, label='f_recon(x)')
plt.legend()
plt.show()