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?
from Example python nfft fourier transform - Issues with signal reconstruction normalization
No comments:
Post a Comment