Skip to content

Improve FFT interpolation (1/2)#116

Open
unalmis wants to merge 10 commits intof0uriest:mainfrom
unalmis:main
Open

Improve FFT interpolation (1/2)#116
unalmis wants to merge 10 commits intof0uriest:mainfrom
unalmis:main

Conversation

@unalmis
Copy link

@unalmis unalmis commented Aug 16, 2025

This is the first of two pull requests to improve the FFT interpolation.

4 goals are achieved by this PR.

  • The real FFT is now used where possible. This is especially useful for 2D transforms. Resolves Use real fft in interpolation #75 .
  • Double the width of the Fourier spectrum is now preserved when interpolating to a less dense grid, at no additional cost. Previously, when the number of requested output samples $n_{\text{out}}$ was less than the number of input samples $n_{\text{in}}$, the interpolation would neglect the $n_{\text{in}} - n_{\text{out}}$ highest frequency spectral coefficients, resulting in lossy interpolation. Now the interpolation is lossless for all $n_{\text{out}} >= \lfloor n_{\text{in}} / 2 \rfloor + 1$. Likewise, if $n_{\text{out}} < \lfloor n_{\text{in}} / 2 \rfloor + 1$, then only the $\lfloor n_{\text{in}} / 2 \rfloor + 1 - n_{\text{out}}$ highest frequency spectral coefficients will be neglected.
  • In the 2D upsampling case, the second transform is now padded only after computing the first transform. This reduces the size of the problem, so the computation is less expensive.
  • In the 2D downsampling case, the second transform is now truncated prior to computing the first transform. This reduces the size of the problem, so the computation is less expensive.

…ncy spectrum

with negligible additional computation.
…transform

is now applied after (before) the first transform is completed.
This reduces the size of the problem, so the computation is less expensive.
return jnp.fft.irfft(c, n, axis=0, norm="forward")

if n < c.shape[0]:
c = c[:n]
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

doesn't this truncation mean that we're still interpolating the band limited signal, not the full one?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

c.shape[0] has size f.shape[0]//2 + 1. so our band limited signal can preserve double the width of the fourier spectrum than before

Comment on lines +57 to +58
x = jnp.linspace(0, 2 * jnp.pi, n, endpoint=False)
x = jnp.exp(1j * (c.shape[0] // 2) * x).reshape(n, *((1,) * (c.ndim - 1)))
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you explain this part? Isn't this just multiplying the output by exp(i*n*x)?

Copy link
Author

@unalmis unalmis Sep 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we transform sum_n^N c_n e^(i n x) for n>=0 to e^iNx//2 sum_k c_n e^(ikx) where k now includes the relevant negative values (that result from pulling the highest frequency outside the sum)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@f0uriest
Copy link
Owner

f0uriest commented Feb 2, 2026

What's the status of this? It looks like there are some merge conflicts. Previously the method was only correct for real valued data, but I've since fixed it so that it works for complex data as well, but I'd be happy to include this as a special case when the data is purely real

@unalmis
Copy link
Author

unalmis commented Feb 2, 2026

but I've since fixed it so that it works for complex data as well, but I'd be happy to include this as a special case when the data is purely real

The improvements in points 3 and 4 in the PR header would still apply to complex transforms.

It looks like there are some merge conflicts.

The merge conflicts are resolved by updating the names of the functions and tests here now that the original ones support complex data. interpax-fft here simply includes #117 with labels rfft instead of fft. As I suspect you may have preferences on API/names, e.g. whether to detect if data is real and call the private funs, I think you should just update and push changes here to your preference.

@unalmis
Copy link
Author

unalmis commented Feb 2, 2026

yea so just add interpax-fft as a depedency and call it if input data is not jnp.complexobj

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

Successfully merging this pull request may close these issues.

Use real fft in interpolation

2 participants