# Sinc Interpolation based on FFT
**URL:** https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512
**Category:** Signal and Image Processing
**Tags:** fftw, interpolations, dsp
**Created:** [December 28, 2020, 1:56pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512 "2020-12-28T13:56:28Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 1:56pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/1 "2020-12-28T13:56:28Z")
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:
```julia
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:
```julia
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:
```julia
#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:
```julia
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
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [December 28, 2020, 2:15pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/2 "2020-12-28T14:15:23Z")
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
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 2:25pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/3 "2020-12-28T14:25:23Z")
Hm, I tried that. First I redefined `xs_low` to be over the full interval:
```julia
xs_low = LinRange(0, 16π, 140)
```
Re-running the interpolation with your hints, returns the following:
```julia
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.
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [December 28, 2020, 2:59pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/4 "2020-12-28T14:59:41Z")
The way I’ve used Fourier interpolation in the past is
```julia
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.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 3:45pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/5 "2020-12-28T15:45:31Z")
Thanks, your suggestion works.
Just for completeness, the full code of my example:
```julia
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!
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 3:48pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/6 "2020-12-28T15:48:09Z")
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].
---
### Author:  [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)
#### Post date: [December 28, 2020, 3:58pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/7 "2020-12-28T15:58:19Z")
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.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 4:00pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/8 "2020-12-28T16:00:30Z")
> [@Arrigo\_Benedetti](#):
>
> 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.
Thanks for the papers, I’ll check them out!
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [December 28, 2020, 4:55pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/9 "2020-12-28T16:55:59Z")
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.
---
### Author:  [@mkitti](https://discourse.julialang.org/u/mkitti)
#### Post date: [December 28, 2020, 6:43pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/10 "2020-12-28T18:43:52Z")
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](https://github.com/JuliaApproximation/ApproxFun.jl)
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [December 28, 2020, 8:06pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/11 "2020-12-28T20:06:25Z")
@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.
```julia
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")
```
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 28, 2020, 8:08pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/12 "2020-12-28T20:08:39Z")
Sorry, you can find it here, it calls a few extra functions:
[https://github.com/JuliaPhysics/PhysicalOptics.jl/blob/edc72ab856aeb3e8cb6282535b72b61a028e6f4d/src/utils.jl#L145](https://github.com/JuliaPhysics/PhysicalOptics.jl/blob/edc72ab856aeb3e8cb6282535b72b61a028e6f4d/src/utils.jl#L145)
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [December 29, 2020, 1:23am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/13 "2020-12-29T01:23:08Z")
@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:
```julia
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
```
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 29, 2020, 11:48am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/14 "2020-12-29T11:48:52Z")
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.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [December 29, 2020, 11:53am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/15 "2020-12-29T11:53:43Z")
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.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 29, 2020, 11:55am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/16 "2020-12-29T11:55:49Z")
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 😂
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 29, 2020, 3:40pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/17 "2020-12-29T15:40:04Z")
I uploaded the code to [GitHub](https://github.com/roflmaostc/FFTInterpolations.jl)
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.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [December 29, 2020, 6:15pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/18 "2020-12-29T18:15:58Z")
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:
```matlab
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.
---
### Author:  [@anon92994695](https://discourse.julialang.org/u/anon92994695)
#### Post date: [December 29, 2020, 7:06pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/19 "2020-12-29T19:06:26Z")
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](https://github.com/caseykneale/ChemometricsTools.jl/blob/397b91b29683700bac101f504773ea967c11f56f/src/Preprocess.jl#L332)
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [December 29, 2020, 7:08pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/20 "2020-12-29T19:08:18Z")
It’s funny that you share this because what my code currently does, is the function `FourierUpsample` below 😃
[Next page](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512.md?page=2)