In short: To process non-uniformly spaced data, apply a Non-Uniform FFT (NUFFT) directly to the irregular samples for frequency-domain analysis, or regrid the samples onto a uniform grid with cubic interpolation and then apply time-domain sinc interpolation (resampleSINC) for a band-limited reconstruction. The article's Python example uses pynufft, scipy interp1d and a sinc kernel on a 256-point two-tone test signal.
When processing non-uniformly spaced data in fields like MRI, radar, or geophysics, it’s crucial to use methods that can handle irregular sampling. In this blog, we discussed two important techniques: Non-Uniform FFT (NUFFT) and Time-domain Sinc Interpolation after regridding the data onto a uniform grid.
Here’s how we integrate both approaches using Python for practical applications.
Step-by-Step Workflow: NUFFT and ResampleSINC Interpolation
- Simulate Non-Uniformly Sampled Data: First, create a sinusoidal signal on a uniform grid, then sample it non-uniformly (like in MRI k-space).
- Apply NUFFT: Use NUFFT to handle the non-uniform data directly.
- Regrid Non-Uniform Data to Uniform Grid: Use interpolation methods like
interp1dto map the data to a uniform grid. - Apply ResampleSINC: Use the
resampleSINCfunction for precise interpolation onto the uniform grid.
Let’s look at the complete Python code:
import numpy as np
import matplotlib.pyplot as plt
from pynufft import NUFFT
from scipy import interpolate
# Step 1: Generate a signal on a uniform grid
N = 256 # Number of points in the signal
x = np.linspace(0, 2*np.pi, N)
signal = np.sin(5*x) + np.sin(15*x) # A test signal with two frequencies
# Step 2: Simulate non-uniform sampling
non_uniform_grid = np.sort(np.random.rand(N) * 2 * np.pi)
# Step 3: Perform NUFFT
nufft = NUFFT()
nufft.plan(non_uniform_grid, (N,), (256,))
signal_non_uniform = np.interp(non_uniform_grid, x, signal)
nufft_transform = nufft.forward(signal_non_uniform)
# Step 4: Regrid non-uniform data onto uniform grid using cubic interpolation
uniform_grid = np.linspace(0, 2*np.pi, N)
regridded_signal = interpolate.interp1d(non_uniform_grid, signal_non_uniform, kind='cubic')(uniform_grid)
# Step 5: Apply ResampleSINC for precise interpolation
def resampleSINC(signal, x, u):
"""
Resample a signal using sinc interpolation in time domain
Parameters:
- signal : original signal on non-uniform grid
- x : non-uniform grid
- u : uniform grid to resample onto
Returns:
- interpolated_signal : the resampled signal on the uniform grid
"""
interpolated_signal = np.zeros_like(u)
for i in range(len(u)):
interpolated_signal[i] = np.sum(signal * np.sinc((u[i] - x) / (x[1] - x[0])))
return interpolated_signal
# Apply ResampleSINC on regridded signal
sinc_interpolated_signal = resampleSINC(regridded_signal, uniform_grid, uniform_grid)
# Step 6: Plot the results
plt.figure(figsize=(12, 6))
# Original and Non-Uniform Sampled Signal
plt.subplot(1, 2, 1)
plt.plot(x, signal, label='Original Signal')
plt.scatter(non_uniform_grid, signal_non_uniform, color='red', label='Non-uniform Samples')
plt.legend()
plt.title('Original and Non-Uniformly Sampled Signal')
# Sinc Interpolated Signal
plt.subplot(1, 2, 2)
plt.plot(x, signal, label='Original Signal')
plt.plot(uniform_grid, sinc_interpolated_signal, label='Sinc Interpolated Signal', linestyle='dashed')
plt.legend()
plt.title('Sinc Interpolation after Regridding')
plt.show()Workflow Breakdown
- Simulate Non-Uniform Sampling: We create a simple sinusoidal signal on a uniform grid and simulate non-uniform sampling, common in MRI k-space.
- Regrid to Uniform Grid: Using cubic interpolation (
interp1dfromscipy), we map the non-uniform samples onto a uniform grid to facilitate subsequent operations like FFT or interpolation. - ResampleSINC: The
resampleSINCfunction is applied after regridding to precisely interpolate the signal onto a uniform grid. This technique uses sinc interpolation, which is ideal for band-limited signals when applied to uniformly spaced data. - Plot Results: We visualize the original, non-uniform samples, and the signal reconstructed using sinc interpolation. As seen in the plot, the sinc-interpolated signal closely follows the original signal, providing an accurate reconstruction.
Conclusion
When working with non-uniformly spaced data, especially in fields like MRI or geophysics, it’s essential to use the right tools to process and reconstruct the signal. This blog demonstrated two key techniques:
- NUFFT allows direct handling of non-uniform data without interpolation, providing accurate frequency-domain transforms.
- ResampleSINC is ideal for signal reconstruction after regridding non-uniform data to a uniform grid, ensuring minimal distortion.
These methods, when used together, offer a powerful solution for handling non-uniformly sampled signals. By first regridding and then applying sinc interpolation, you ensure your final signal is both smooth and faithful to the original data.
Reference
Dr. Erol Kalkan, P.E. (2024). Time-domain Sinc Interpolation (Resampling) (https://www.mathworks.com/matlabcentral/fileexchange/59027-time-domain-sinc-interpolation-resampling), MATLAB Central File Exchange. Retrieved September 17, 2024.
If you’re dealing with non-uniform sampling in your field, contact QuakeLogic today for advanced solutions at info@quakelogic.net.
Let us show you how our state-of-the-art technologies can transform your signal processing workflows.
Key Takeaways
- NUFFT (pynufft) transforms non-uniformly sampled data to the frequency domain without first interpolating.
- scipy.interpolate.interp1d with kind='cubic' regrids irregular samples onto a 256-point uniform grid.
- resampleSINC sums signal * sinc((u – x) / dx) over all samples, an exact reconstruction for band-limited uniformly sampled data.
- The example signal is sin(5x) + sin(15x) on 0 to 2π with randomly placed non-uniform sample points.
- The method follows Erol Kalkan's Time-domain Sinc Interpolation (Resampling) submission on MATLAB File Exchange.
Frequently Asked Questions
A Non-Uniform FFT (NUFFT) computes the Fourier transform of data whose sample points are irregularly spaced, which a standard FFT cannot accept. Use NUFFT when the data arrives on a non-uniform grid, as in MRI k-space, radar or geophysical surveys, and you need a frequency-domain result without first interpolating. In Python the pynufft package plans the transform from the sample positions and grid size, then calls forward() on the samples.
Sinc interpolation reconstructs a band-limited signal exactly from uniformly spaced samples by summing each sample multiplied by a sinc kernel centred on it. Because the kernel assumes uniform spacing, the workflow first regrids the irregular samples with cubic interpolation, then applies the resampleSINC function, which computes sum(signal * sinc((u – x) / dx)) for every output point u. The result follows the original signal closely with minimal distortion.
Time-domain sinc interpolation is only exact when the input samples lie on a uniform grid, so the irregular samples must be mapped onto one first. scipy.interpolate.interp1d with kind='cubic' provides a smooth, fast approximation of the signal at the uniform grid points; resampleSINC then refines that estimate with a band-limited kernel, giving a reconstruction that is both smooth and faithful to the original data.
The example uses NumPy for arrays and the sinc function, Matplotlib for plotting, pynufft for the NUFFT class with its plan and forward calls, and SciPy's interpolate module for cubic regridding. The resampleSINC function itself is a short pure-NumPy loop, so no additional dependency is required for the sinc stage of the workflow.
Yes. Geophysics is one of the fields the method targets, alongside MRI and radar. Field sensors with irregular timestamps, gaps or uneven station spacing produce non-uniform samples that break standard FFT-based processing. Regridding followed by sinc resampling restores a uniform time series suitable for spectral analysis, while NUFFT gives a direct spectrum when interpolation should be avoided.
Related Articles
Need help with a project? QuakeLogic designs and supplies seismic monitoring, earthquake early warning, structural health monitoring, infrasound, and shake table systems. Browse QuakeLogic products or contact our engineering team with your site, structure, or test requirements.
