Sinc Interpolation based on FFT

roflmaostc · December 28, 2020, 1:56pm

Hey,

I was trying to implement a sinc interpolation based on FFTs.

One way (as far as I know) is to zero pad the frequency spectrum and calculate the IFFT of the spectrum to get a sinc interpolated signal.

I ended up with the following implementation:

function sinc_interpolate(arr::AbstractArray{T}, new_size) where T<:Complex
    arr_f = fftshift(fft(arr))
    out_f = zeros(eltype(arr_f), new_size)
    center_set!(out_f, arr_f)
    out = ifft(ifftshift(out_f)) ./ length(arr) .* length(out_f)
    return real(out) 
end

center_set! is a function which sets the signal to the center of the array (in the FFT convention).
To compare the method, I also implemented a slow, sum based sinc interpolation:

function sinc_interpolate_sum(arr, new_size)
    out = zeros(eltype(arr), new_size)
    T = 1 / size(arr)[1]
    t_arr = LinRange(0, 1.0 - T, size(out)[1])
    for (j, t) in enumerate(t_arr)
        v = zero(eltype(arr))
        for n = 1:size(arr)[1]
            v += arr[n] * sinc((t - (n-1) * T) / T)
        end
        out[j] = v
    end
    return out
end

The FFT based method basically works and I’m familiar with the fact, that this kind of interpolation assumes a infinite signal and a finite frequency support.

However, I observe certain artifacts with simple examples:

#prepare x pos and function
xs_low = LinRange(0, 16π - 16π/140, 140)
f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)

#calculate low resolution
arr_low = f.(xs_low)

#interpolation
N = 1000
xs_interp = LinRange(0, xs_low[end], N)
arr_interp = sinc_interpolate(arr_low, N)
arr_interp_s = sinc_interpolate_sum(arr_low, N)

We can plot everything:

scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="data")
plot!(xs_interp, arr_interp, label="FFT sinc interp")
plot!(xs_interp, arr_interp_s, label="sum sinc interp")

We see that the FFT based curve (orange) seems to have a slight, growing offset despite I’ve chosen the data arr_low to be periodic.

My questions would be:

Is this shift due to the finite signal?
Can we fix this shift somehow? For example, by stretching the signal, it fits better.
Is this even a proper sinc interpolation? Of course, I was searching in the web, but I often find vague, unprecise information People mention sometimes the aliased sinc or dirichlet kernel. Can anybody refer a good theoretical resource?

Thanks,

Felix

antoine-levitt · December 28, 2020, 2:15pm

Haven’t looked at the code closely, but usually when working with periodic functions you don’t want to include both endpoints in the interval

roflmaostc · December 28, 2020, 2:25pm

Hm, I tried that. First I redefined xs_low to be over the full interval:

xs_low = LinRange(0, 16π, 140)

Re-running the interpolation with your hints, returns the following:

xs_interp = LinRange(xs_low[2], xs_low[end-1], N)
arr_interp = sinc_interpolate(arr_low[2:end-1], N)

So there is still this offset.

antoine-levitt · December 28, 2020, 2:59pm

The way I’ve used Fourier interpolation in the past is

x_low = range(a, b, length=N_low+1)[1:N_low]
x_high = range(a, b, length=N_high+1)[1:N_high]
f_low = f.(x_low)
# fft f_low, pad with zeros to N_high, ifft back, get f on `x_high`

The math of this is that, by unitarity, the fft of f_low gets you the Fourier coefficients of the only Fourier series that passes through the points (x_low, f_low). Then padding with zeros and iffting back is simply the evaluation of that Fourier series. In order for the formulas to match, the x_low and x_high have to be as above.

roflmaostc · December 28, 2020, 3:45pm

Thanks, your suggestion works.

Just for completeness, the full code of my example:

function sinc_interpolate_sum(arr, new_size)
    out = zeros(eltype(arr), new_size)
    T = 1 / size(arr)[1]
    t_arr = range(0, 1.0, length=size(out)[1]+1)[1:size(out)[1]]
    for (j, t) in enumerate(t_arr)
        v = zero(eltype(arr))
        for n = 1:size(arr)[1]
            v += arr[n] * sinc((t - (n-1) * T) / T)
        end
        out[j] = v
    end
    return out
end

begin
	N_low = 128
	x_min = 0.0
	x_max = 16π
	
	xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
	f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
	arr_low = f.(xs_low)
end

begin
	N = 1000
	xs_interp = range(x_min, x_max, length=N+1)[1:N]
	arr_interp = sinc_interpolate(arr_low, N)

	N2 = 1000
	xs_interp_s = range(x_min, x_max, length=N2+1)[1:N2]
	arr_interp_s = FFTInterpolations.sinc_interpolate_sum(arr_low, N2)
end

begin
	scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="data")
	plot!(xs_interp, arr_interp, label="FFT sinc interp")
	plot!(xs_interp_s, arr_interp_s, label="sum sinc interp")	
end

Thanks for your help!

roflmaostc · December 28, 2020, 3:48pm

Maybe a extra question you know:

Sometimes I find the information, that for a even FFT, people not only zero pad, but also modify the first non-zero positive frequency to achieve the Hermitian property.
If the padded spectrum would be: [0 0 C* B* A B 0 0] they instead calculate the IFT of [0 0 C*/2 B* A B C*/2 0].

Arrigo_Benedetti · December 28, 2020, 3:58pm

When doing sinc interpolation with the FFT there are some subtle technicalities that have to be considered. I don’t remember the math off the top of my head, it has to do with a re-scaling of the first and last samples of the sequence in the time domain. All the details are in these paper:

Schanze, Thomas. “Sinc interpolation of discrete periodic signals.” IEEE Transactions on Signal Processing 43.6 (1995): 1502-1503.

Candocia, Frank, and Jose C. Principe. “Comments on” Sinc interpolation of discrete periodic signals." IEEE Transactions on signal processing 46.7 (1998): 2044-2047.

roflmaostc · December 28, 2020, 4:00pm

Thanks for the papers, I’ll check them out!

antoine-levitt · December 28, 2020, 4:55pm

It’s equivalent to taking the real part of the output array, so there’s no specific need for it unless you use the real versions of the fft. Should only be necessary if super sampling from an odd to an even grid or something like that.

mkitti · December 28, 2020, 6:43pm

You may want to try this on a windowed signal where the signal at both ends is attenuated to zero by some smooth function. The easiest would be a Gaussian window. The periodization of the windowed signal should provide better results.

For a finite bounded signal like this, something like Chebyshev polynomials would provide a better basis than periodic functions. For example see:

https://github.com/JuliaApproximation/ApproxFun.jl

rafael.guerra · December 28, 2020, 8:06pm

@roflmaostc, you did not seem to provide the center_set! function to reproduce the results.
In code snippet below the sinc interpolation function uses fftshift. It seems to work fine for even number of samples (both input and output). Did not figure out now how to extend this function to arbitrary parity.

using FFTW, Plots

function sinc_interpolate(arr, new_size::Int)   # for IO with even number of samples
    n = length(arr)
    arr_f = fft(arr)
    arr_f[n÷2+1] = 0.5*(arr_f[n÷2] + arr_f[n÷2+2])  #to minimize FFT artifacts
    arr_f = fftshift(arr_f)
    k = (new_size - n + 1) ÷ 2
    arr_f = [zeros(k,1); arr_f; zeros(k,1)]
    return new_size/n * ifft(fftshift(arr_f))[1:new_size] 
end

# Signal at low sampling rate
n = 128  # even number of samples
t0, t1 = 0, 16π
xs_low = LinRange(t0, t1, n+1)[1:n]
f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
arr_low = f.(xs_low)  # evaluate f at low sampling rate

# Sinc interpolation
new_size = 500   # even number of samples
xs_interp = LinRange(t0, t1, new_size + 1)[1:new_size]
arr_interp = real(sinc_interpolate(arr_low, new_size))

scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="input data")
plot!(xs_interp, real.(arr_interp), label="FFT sinc interp")

sinc_interpolation

roflmaostc · December 28, 2020, 8:08pm

Sorry, you can find it here, it calls a few extra functions:
https://github.com/JuliaPhysics/PhysicalOptics.jl/blob/edc72ab856aeb3e8cb6282535b72b61a028e6f4d/src/utils.jl#L145

rafael.guerra · December 29, 2020, 1:23am

@roflmaostc, thank you for sharing the functions to get the indices around the center for general case. Based on your code, and for 1D arrays, a standalone sinc_interpol() function may be written as follows:

function sinc_interpol(arr, ns::Int)
    n = length(arr)
    arr_f = fftshift(fft(arr))
    c = ns + isodd(n) - isodd(ns)
    x = (c - n) ÷ 2
    a, b = 1 + x, c - x
    out = fill(0.0im, ns, 1)
    out[a:b] = arr_f
    return ifft(ifftshift(out)) * ns/n
end

roflmaostc · December 29, 2020, 11:48am

A short remark, ifftshift is always needed not only for ns being odd.

I’m probably going to push the working code for N-d arrays soon to a repository. But before I want to check the referenced papers above.

rafael.guerra · December 29, 2020, 11:53am

Sorry, I agree that the remark was not clear, what I meant is why ifftshit is needed and not fftshift. From the help file: “For even N this is equivalent to swapping the first and second halves, so fftshift and ifftshift are the same.”. PS: I have edited it now.

roflmaostc · December 29, 2020, 11:55am

I guess it’s always better to use fftshift to shift the center frequency to the center position.
And ifftshift to shift the center frequency back from the center to the top left corner.
Once memorized, you’re never confused about that anymore

roflmaostc · December 29, 2020, 3:40pm

I uploaded the code to GitHub

Maybe we can join efforts, the basic code already works, but so far I haven’t implement the different sum algorithms from the paper above.

The examples folder contains the plot example from above.

RoyiAvital · December 29, 2020, 6:15pm

If I get you correctly, what you want is to interpolate the signal in the spectrum.

Many years ago I derived it. I don’t have the code in Julia, Back then I only used MATLAB.
But it will be easy to convert:

function [ vX ] = SincInterpolationDft( vXDft, numSamplesTime )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = SincInterpolationDft( vXDft, numSamplesTime )
% Applies Sinc interpolation to the input signal while preserving the
% Parseval Theorem energy.
% Input:
%   - vXDft             -   DFT Samples.
%                           The DFT samples of the vector to be
%                           interpolated.
%                           Structure: Vector (numSamples x 1).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - numSamplesTime    -   Number of Samples in Time.
%                           The number of samples of the output in time.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
% Output:
%   - vX                -   Output Vector.
%                           The interpolated vector samples in time.
%                           Equivalent to be interpolated by Sinc (Actually
%                           by Dirichlet Kernel).
%                           Structure: Vector (numSamplesTime x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  A.
% Remarks:
%   1.  The samples are assumed to be in [0, 2 pi) range (Namely the output
%       of `fft()` with no `fftshift()`).
%   2.  The function takes care of the even and odd case of the samples. It
%       assumes the original signal was real hence keeps on hermitian
%       symmetry and the Parseval Theorem.
%   3.  In order to utilize the interpolation the number of the samples in
%       time must be larger than the input. Namely 'numSamplesTime >
%       size(vXDft, 1)'.
% TODO:
%   1.  
%   Release Notes:
%   -   1.0.001     29/12/2020  Royi Avital
%       *   Updated to modern MATLAB.
%   -   1.0.000     20/01/2009  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    vXDft (:, 1) {mustBeNumeric}
    numSamplesTime (1, 1) {mustBeNumeric, mustBeReal, mustBePositive, mustBeLonger(numSamplesTime, vXDft)}
    % numSamplesTime (1, 1) {mustBeNumeric, mustBeReal, mustBePositive}
end

numSamples = length(vXDft);

if(numSamplesTime <= numSamples)
    vX = ifft(vXDft, numSamplesTime);
    return;
end

if (mod(numSamples, 2) ~= 0) % Odd number of samples
    vXDftInt = (numSamplesTime / numSamples) * [vXDft(1:ceil(numSamples / 2)); zeros(numSamplesTime - numSamples, 1, 'like', vXDft); vXDft((ceil(numSamples / 2) + 1):numSamples)];
else % Even number of samples -> Special Case
    vXDftInt = (numSamplesTime / numSamples) * [vXDft(1:(numSamples/2)); vXDft((numSamples / 2) + 1) / 2; zeros(numSamplesTime - numSamples - 1, 1, 'like', vXDft); vXDft((numSamples / 2) + 1) / 2; vXDft(((numSamples / 2) + 2):numSamples)];
end
  
vX = ifft(vXDftInt);


end


function [ ] = mustBeLonger( numSamples, vX )
    if(numSamples < size(vX, 1))
        eid = 'mustBeLonger:numSamplesMustBeLonger';
        msg = 'The number of smples of the output must be at lest the number of samples of the input';
        throwAsCaller(MException(eid, msg));
    end
    
    
end

The idea is to keep the conjugate symmetry and the Parseval Theorem.
Both means that in the case of even number of samples you split the element of the highest frequency.
If you impose both constraints, you’d see this is the solution.

Somewhere I also have the case for less samples. I have to find it.

anon92994695 · December 29, 2020, 7:06pm

Probably not what you want - but I remember stumbling on this a while ago, cool 1 liner

https://github.com/caseykneale/ChemometricsTools.jl/blob/397b91b29683700bac101f504773ea967c11f56f/src/Preprocess.jl#L332

roflmaostc · December 29, 2020, 7:08pm

It’s funny that you share this because what my code currently does, is the function FourierUpsample below