Handling Non-Uniformly Spaced Data Using NUFFT and Sinc Interpolation

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

1. Simulate Non-Uniformly Sampled Data: First, create a sinusoidal signal on a uniform grid, then sample it non-uniformly (like in MRI k-space).
2. Apply NUFFT: Use NUFFT to handle the non-uniform data directly.
3. Regrid Non-Uniform Data to Uniform Grid: Use interpolation methods like interp1d to map the data to a uniform grid.
4. Apply ResampleSINC: Use the resampleSINC function 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

1. Simulate Non-Uniform Sampling: We create a simple sinusoidal signal on a uniform grid and simulate non-uniform sampling, common in MRI k-space.
2. Regrid to Uniform Grid: Using cubic interpolation (interp1d from scipy), we map the non-uniform samples onto a uniform grid to facilitate subsequent operations like FFT or interpolation.
3. ResampleSINC: The resampleSINC function 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.
4. 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

Related Articles

• Calculating Background Seismic Noise Levels: A Step-by-Step Guide
• Linearity, Repeatability and Phase Lag in Digital Sensors
• How to Compute the PGA, PGV, PGD and Response Spectrum of Ground Motion Using Python
• The Doppler Effect: A Powerful Tool for Structural Health Monitoring

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.