Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Equispaced data with gaps vs Non-equispaced data #20

Open
wgprojects opened this issue Aug 30, 2022 · 7 comments
Open

Equispaced data with gaps vs Non-equispaced data #20

wgprojects opened this issue Aug 30, 2022 · 7 comments

Comments

@wgprojects
Copy link

Can anyone help clear up my misunderstanding?

I used the forward transform to convert a sine wave to frequency domain, then adjoint transform to go back to time domain.
If the data is equispaced, this works well (not pictured).

If the data is equispaced but with even a small gap (<1% missing) the reconstruction is smooth, but with incorrect amplitude:
Gap

If the data is non-equispaced, the reconstruction is very poor. The Fourier coefficients appear totally random - not an ideal spike, and not spectral leakage.
NEq

Is this technique not intended for just these applications? I believe it is and that I've made a mistake.

`import nfft
import numpy as np
import matplotlib.pyplot as plt

freq = 13 #Rad/s
N = 300

x = np.linspace(-0.5, +0.5, N)
f = np.sin(freq*x)

#Equally spaced data with one gap:

gapsize = 2
x=np.delete(x, slice(N//2,N//2+gapsize), axis=0)
f=np.delete(f, slice(N//2,N//2+gapsize), axis=0)

N = len(x)
k = -N//2 + np.arange(N)

f_k = nfft.ndft(x, f)

plt.figure(2)
plt.scatter(k, f_k.real, label="Real")
plt.scatter(k, f_k.imag, label="Imag")
plt.legend()

f_hat = nfft.ndft_adjoint(x, f_k, N)/N
plt.figure(1)
plt.scatter(x, f, s=20, marker='o', label="Data with gaps")
plt.scatter(x, f_hat, s=35, marker='v', label="NDFT fit")

plt.legend()

#Randomly spaced data
x = 0.5 + np.random.rand(N)
f = np.sin(freq*x)

N = len(x)
k = -N//2 + np.arange(N)

f_k = nfft.ndft(x, f)

plt.figure(4)
plt.scatter(k, f_k.real, label="Real")
plt.scatter(k, f_k.imag, label="Imag")
plt.legend()

f_hat = nfft.ndft_adjoint(x, f_k, N)/N
plt.figure(3)
plt.scatter(x, f, s=20, marker='o', label="Data non-equispaced")
plt.scatter(x, f_hat, s=35, marker='v', label="NDFT fit")

plt.legend()
plt.show()
`

@PNeigel
Copy link

PNeigel commented Jan 3, 2023

It seems like ndft_adjoint is actually the function to go from signal into frequency space, and ndft to go from frequency into signal space.

freq = 3
N_samples = 500

x = np.linspace(0, 1, N_samples) # 1 is arbitrary here. This is the time range of your data
f = np.sin(freq*2*np.pi*x)

# Equally spaced data

plt.scatter(x,f)
plt.title("Equally spaced signal")
plt.show()

N_freq = 100
k = -N_freq//2 + np.arange(N_freq)

f_k = nfft.ndft_adjoint(x, f, N_freq)

plt.figure(2)
plt.plot(k, np.abs(f_k), label="Absolute")
plt.legend()
plt.title("Frequency spectrum for equally spaced signal")
plt.show()

print("Maximum at ", k[np.argmax(np.abs(f_k))])

f_hat = nfft.ndft(x, f_k)/N_samples
plt.figure(1)
plt.scatter(x, f, s=35, marker='o', label="Data with gaps")
plt.scatter(x, f_hat, s=15, marker='v', label="NDFT fit")
plt.legend()
plt.title("Reconstruction for equally spaced signal")
plt.show()

#Equally spaced data with one gap:

gapsize = 30
x=np.delete(x, slice(N_samples//2,N_samples//2+gapsize), axis=0)
f=np.delete(f, slice(N_samples//2,N_samples//2+gapsize), axis=0)

N_freq = 100
k = -N_freq//2 + np.arange(N_freq)

f_k = nfft.ndft_adjoint(x, f, N_freq)

plt.figure(2)
plt.plot(k, np.abs(f_k), label="Absolute")
plt.legend()
plt.title("Frequency spectrum for equally spaced signal with gap")
plt.show()

print("Maximum at ", k[np.argmax(np.abs(f_k))])

f_hat = nfft.ndft(x, f_k)/(N_samples-gapsize)
plt.figure(1)
plt.scatter(x, f, s=35, marker='o', label="Data with gaps")
plt.scatter(x, f_hat, s=15, marker='v', label="NDFT fit")
plt.legend()
plt.title("Reconstruction for equally spaced signal with gap")
plt.show()

#Randomly spaced data
x = np.random.rand(N_samples)
f = np.sin(freq*2*np.pi*x)

N = len(x)
k = -N_freq//2 + np.arange(N_freq)

f_k = nfft.ndft_adjoint(x, f, N_freq)

plt.figure(4)
plt.plot(k, np.abs(f_k), label="Absolute")
plt.title("Frequency spectrum for randomly spaced signal")
plt.legend()
plt.show()

print("Maximum at ", k[np.argmax(np.abs(f_k))])

f_hat = nfft.ndft(x, f_k)/N_samples
plt.figure(3)
plt.scatter(x, f, s=20, marker='o', label="Data non-equispaced")
plt.scatter(x, f_hat, s=35, marker='v', label="NDFT fit")
plt.title("Reconstruction for randomly spaced signal")
plt.legend()
plt.show()

result:

grafik
grafik
grafik
grafik

@evileleven
Copy link

I am also quite confused about whether nfft_adjoint or nfft correspond to the Fourier transform? which is the inverse Fourier transform.

@jakevdp
Copy link
Owner

jakevdp commented Oct 15, 2024

Depending on what Fourier transform convention you use, the equivalent might be nfft_adjoint or nfft.

@evileleven
Copy link

@jakevdp thank you for your reply!
If I want to do nonuniform DFT (the sample interval for time data is change),where x[n] is the time data,X[k] is the frequency data.
image
Should I use nfft_adjoint or nfft?
I think I should use nfft, but the result seems tell me to use nfft_adjoint.

@jakevdp
Copy link
Owner

jakevdp commented Oct 29, 2024

The formula you quote looks like a standard DFT (both n and k are on a regular grid) so I would suggest a standard FFT for that case.

If you have non-uniform data or frequencies, you can use the nfft in this package, and the mathematical form of the forward and adjoint transforms can be seen here: https://github.com/jakevdp/nfft?tab=readme-ov-file#about. You should use the form that most closely matches the quantity you wish to compute.

You mention that "the sample interval for time data is change", and I interpret that to mean that you want the output to be on a regular frequency grid, so that suggests the adjoint transform is most appropriate. Additionally, you'll have to correct the output for the difference in phase convention used in the DFT formula you mention.

Hope that helps!

@wgprojects
Copy link
Author

wgprojects commented Oct 29, 2024 via email

@jakevdp
Copy link
Owner

jakevdp commented Oct 29, 2024

This package provides a fast implementation of two particular summations. Honestly, if the formulas for those summations don’t make sense to you, this might not be the right tool for you to use.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants