Sinc Interpolation based on FFT

I would argue that the correct symmetry handling is exactly that which results in zeroing out the imaginary part of the result, without changing its real part. The geometrical way of seeing this is that you’ve got the subspace X of purely real real-space signals, which is dual by the FFT to the subspace Y of frequency-space sequences satisfying a particular symmetry. If you’ve got a frequency-space signal y that doesn’t belong to Y, you can project it in that subpsace, ie find the ytilde in Y that is closest to y in the least-squares. By Parseval, that operation corresponds in real space to projecting in the least-squares sense the ifft of y on X, which is exactly taking the real part. Again, the nice thing with this is that 1) you don’t have to make a choice for how you fix the nyquist components, that choice is given by the least-squares projection 2) you don’t even have to make explicit what Y is, or what the projection on Y is. Of course you could write it down explicitly, but why would you bother doing that?

I’m bothering because taking the real part of the iffting result produces worse results than my current routine:

begin
	using Revise, FFTInterpolations, FFTW
end

begin
	x = range(-10.0, 10.0, length=129)[1:end-1]
	y = x'
	arr = abs.(x) .+ abs.(y) .+ sinc.(sqrt.(x .^2 .+ y .^2))
end

function ds_r(y)
	a = fftshift(fft(y))
	b = a[1:end-1, 1:end-1]
	b[1:end, 1] .*= 2
	b[1, 1:end] .*= 2
	b[1, 1] /= 2
	return real(ifft(ifftshift(b))) .* length(b) ./ length(y)
end
# 19167.48880005109
sum(abs2.(arr))

sum(abs2.(arr .- ds_r(sinc_interpolate(arr, (129, 129)))))
# 1.1928437627612847e-5

sum(abs2.(arr .- downsample(sinc_interpolate(arr, (129, 129)), (128, 128))))
# 5.592278121264315e-27

However, what is more in accordance, if I modify ds_r:

function ds_r(y)
	a = fftshift(fft(y))
	b = a[1:end-1, 1:end-1]
	b[1:end, 1] .*= 2
	b[1, 1:end] .*= 2
	b[1, 1] /= 2
	return real(ifft(ifftshift(b))) .* length(b) ./ length(y)
end

sum(abs2.(arr .- ds_r(sinc_interpolate(arr, (129, 129)))))
# 5.483458710517112e-27

sinc_interpolation makes the array hermitian and divides the boundaries by 2 (as @RoyiAvital did on StackExchange.

OK so your point is that upsampling then downsampling a real signal has to produce the same result? Indeed the naive scheme of ignoring symmetry and then taking real parts does not preserve this property: if eg upsampling from 2 to 3, the FFT [a b] will get upsampled to [a b 0], but then taking the real part will get [a b/2 b*/2], which will get downsampled to [a b/2]. I don’t think it’s a big deal: the error is of the order of the amplitude of the Nyquist component, which should be small - note that in your experiment you interpolated a function with discontinuous derivatives, hence the relatively large error even at high N. If you really really want to preserve that property indeed it makes sense to look for fancy schemes, but I would be surprised if you didn’t lose other properties in the process.

Yeah after some time I thought I would be a nice property. Otherwise it is hard to test the functions.

The discontinuity is unfortunate, but the result was similar for the pure sinc.

I’m satisfied with the current state of FFTInterpolations.jl
And I leave it like this for now.

Thanks for your input!

If you’re downsmapling from N to M<N samples and M is an even number all needed is to keep the first (In the fftshift() grid) M/2 element (At the index -M/2), multiply it by 0.5 and take only its real part. You remove out the other new nyquist sample (At the M/2 index) and all other smaples with indes greater than M/2 then do regular fft() and the result should be real (Up to floating point errors).

Just as in my code, just the real part of the negative nyquist sample.

@Arrigo_Benedetti, thanks for sharing the Candocia and Principe (1998) reference. It addresses 1D-sinc interpolation of periodic band-limited complex series but it might be extended to more dimensions. It may complement @roflmaostc’s FFTResampling.jl package, by allowing interpolating a coarsely regularly sampled array at a set of arbitrary points within the support of the input.

Some examples below for 1D-sinc interpolation of real and complex series at arbitrary points:

The example below of 2D-sinc interpolation uses two-passes, along x and y, of the above Candocia’s 1D-sinc interpolation. It includes also a nice resampling result from FFTResampling.jl:

The basic Julia code with examples and annotations takes some 150 lines and for that reason was not posted it here. Feel free to ping me if it might be of interest.

3 Likes

I was also considering to rename it and you gave men the final push :wink:
FFTResampling.jl

2 Likes

Very nice work! I love the 3D complex plot …

3 Likes

A really neat package that I now use to interpolate my periodic functions (not in the field of signal processing though). I’d say it is a perfect fit for my package which also uses FFTW internally. Thank you for your great work!

To my use-case, I will do either interpolation or downsampling in a single optimization run. I notice that sinc_interpolate and downsample methods have similar (or same?) interfaces. I thus wonder is it possible to have a, for example, resample method that combines the two? Currently, I use

resample = N_new > N_old ? sinc_interpolate : downsample
A_new = resample(A_old, N_new) 
3 Likes

Yes, thanks for your suggestion that’s definitely possible. It’s just happened to be so without much considerations :grinning:

I’ll suggest you leave the workaround and I will release a breaking 0.3.0 soon!

4 Likes

I thought about your wish.
We can even generalize resample to handle both down and upsampling at the same time!

julia> x = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> size(x)
(3, 3)

julia> resample([1 2 3; 4 5 6; 7 8 9], (4, 2))
4×2 Array{Float64,2}:
 1.0      3.0
 2.26795  4.26795
 7.0      9.0
 5.73205  7.73205

The new code is here

It is not very elegant, but I basically do first a upsampling of all dimensions which are larger, and then I do the final downsampling step. In that way we can mix the operations in one step.

I’ll try to clean everything (examples etc.) and then tag a new version by the end of the week.

3 Likes