# 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:** 3

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 7, 2021, 9:52pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/41 "2021-01-07T21:52:51Z")

</div>

By “Nyquist Sample” you mean the highest frequency bin?

Yeah, Julia doesn’t have the `symmetric` flag. I guess you are referring to `rfft` instead?

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 7, 2021, 10:00pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/42 "2021-01-07T22:00:39Z")

</div>

Nyquist Sample only exist in the case the number of samples in Frequency Domain is even. Then it is the `n / 2` sample (Before `fftshift()`) or the 1st (After `fftshift()`).  
If you don’t have symmetric option, then just multiply it by 0.5 and take its real part.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 8, 2021, 9:23am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/43 "2021-01-08T09:23:24Z")

</div>

I think the using symmetric (`rfft`) is probably the better choice because by design it does the right thing, doesn’t it?

I was wondering about a complex array. Wouldn’t it be also better to use `rfft` for imaginary part and real part separately and add them afterwards back together?  
Because for a purely real/complex array we know the symmetry property and we know that `rfft` does the right thing.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 8, 2021, 12:24pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/44 "2021-01-08T12:24:25Z")

</div>

Puh, this one never ends.

Using `rfft` and `irfft` also induces some problems because the `irfft` expects that the input is basically hermitian (if you would unfold it to the original size). If you ensure this, you get entirely wrong results in some rows.

In this case the first row is mostly wrong.

```julia
julia> z
5×8 Array{Complex{Float64},2}:
  0.809667-0.223829im -0.146846+0.483231im 0.126273-0.673676im -0.0995014-0.652765im 1.08065+0.179737im -0.0545377-1.29798im 0.291862+0.335662im 0.244291-0.939019im
 0.0796073-0.551711im -0.261923-1.15531im -0.565023-0.699141im -0.374933-1.15981im -0.308387-1.51755im -0.194571-1.77467im -0.912418+0.388777im 0.0659049+0.528615im
  -0.49447+0.298023im -0.905227+0.545565im -1.10754+0.15279im -0.776708+0.0774121im 0.647339-0.294794im 0.426044-0.803606im 0.393879-0.603508im -0.653537+0.483806im
  0.521922-0.491875im -0.585717+1.07584im 1.85847-0.241757im -1.08196+0.141852im 0.990565-0.123956im 1.42857-0.675697im -0.0593938-0.513508im 0.353353+0.538698im
  0.672973+0.311109im -0.134737-0.770579im 0.571894-1.15421im 1.0748-0.245317im -0.226026+0.132937im -0.424684+0.98478im -0.389344-0.941847im 0.866529-0.262219im

julia> rfft(irfft(z, 8))
5×8 Array{Complex{Float64},2}:
  0.809667+0.0im 0.0487227+0.711125im 0.209068-0.504669im -0.0770195+0.322609im 1.08065+0.0im -0.0770195-0.322609im 0.209068+0.504669im 0.0487227-0.711125im
 0.0796073-0.551711im -0.261923-1.15531im -0.565023-0.699141im -0.374933-1.15981im -0.308387-1.51755im -0.194571-1.77467im -0.912418+0.388777im 0.0659049+0.528615im
  -0.49447+0.298023im -0.905227+0.545565im -1.10754+0.15279im -0.776708+0.0774121im 0.647339-0.294794im 0.426044-0.803606im 0.393879-0.603508im -0.653537+0.483806im
  0.521922-0.491875im -0.585717+1.07584im 1.85847-0.241757im -1.08196+0.141852im 0.990565-0.123956im 1.42857-0.675697im -0.0593938-0.513508im 0.353353+0.538698im
  0.672973+0.0im 0.365896-0.25418im 0.0912752-0.10618im 0.32506-0.615049im -0.226026+0.0im 0.32506+0.615049im 0.0912752+0.10618im 0.365896+0.25418im

```

So these functions are not the inverse of each other.

It basically ends at the same point, that you somehow need to make them “hermitian”. I already have this one, but for full fft(arr) arrays. For the `rfft` it seems just be more complex.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 8, 2021, 9:06pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/45 "2021-01-08T21:06:52Z")

</div>

Maybe, I finally reached a working state which I’m satisfied with.

## Sinc Interpolation with FFT

I’m always using `fft` for both real and complex signals. For real signals, in the even length signal case, I artificially change certain frequencies on the positive side, to create an array which has the hermitian symmetry. By doing so, the result from iffting is purely real.

## Sinc Downsampling

I’m always using `fft` for both real and complex signals. In downsampling where my output size is even in an dimension, I artificially change the highest negative frequencies (where the positives partners are missing) to obey the symmetries so that iffting returns a purely real signal.

All in all, upsampling and downsampling are approximately inverse to each other and recover the initial signal pretty well.

I also registered that package and it’ll be merged in the [next days](https://github.com/JuliaRegistries/General/pull/27600)

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [January 8, 2021, 9:28pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/46 "2021-01-08T21:28:46Z")

</div>

I still don’t understand why you bother. Just ignore the symmetries (they are bothersome to deal with, esp in higher dimensions), and if the input to the function was real, just discard the imaginary part after up/down sampling.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 8, 2021, 11:15pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/47 "2021-01-08T23:15:35Z")

</div>

Hm, didn’t feel right to do so.

With the current routine, upsamling and downsampling are better inverse to each other in comparison to the routine where we take always the real part after iffting.

Nevertheless , thanks for your valuable input, I definitely learnt something!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2021, 1:59am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/48 "2021-01-09T01:59:24Z")

</div>

> [@antoine-levitt](#):
>
> Just ignore the symmetries (they are bothersome to deal with, esp in higher dimensions), and if the input to the function was real, just discard the imaginary part after up/down sampling.

If the resulting spurious imaginary parts are large, they probably indicate a bug in your algorithm. It’s pretty easy to get FFT-based interpolation (or differentiation etc.) wrong.

`norm(imag(output)) / norm(real(output))` should be nearly on the order of machine precision.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [January 9, 2021, 7:50am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/49 "2021-01-09T07:50:39Z")

</div>

Not in this very mparticular case (downsampling from odd to even size). Of course the resulting imaginary part should be small (or your downsampling will be bad anyway and you should have filtered it before) but not of order of the machine precision.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 9, 2021, 11:20am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/50 "2021-01-09T11:20:47Z")

</div>

This is the wrong policy in my opinion.  
But the code I supplied solves the case for downsmapling.

I really don’t understand what’s missing.

Not to say you can always do the downsmapling by upsampling and decimation.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 9, 2021, 11:30am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/51 "2021-01-09T11:30:36Z")

</div>

@RoyiAvital It does maybe in 1D, but as far as I tested this, it doesn’t generalize to N dimensions. After iffting, the result contains a non-zero imaginary part. Using a `symmetric` FFT (as you did) provides real results, but at least with FFTW (via `rfft`) the real results are messed up a lot (low and high frequencies in the first dimension).

And downsampling by decimation introduces aliasing because you simply undersample a high frequency signal. Taking a frequency window from a FFT spectrum is not the same as decimation because it is alias free.  
(by decimation you mean taking every N-th sample?)

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 9, 2021, 11:47am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/52 "2021-01-09T11:47:06Z")

</div>

What I mean:  
`rfft` requires a special symmetry to produce a correct real result. If you don’t provide this symmetry, you really mess up low and high frequencies:

The symmetry can be seen here:

```julia
julia> rfft(randn((6,6)))
4×6 Array{Complex{Float64},2}:
  4.36498+0.0im 5.59917+4.04078im 2.42256+1.53944im 0.583178+0.0im 2.42256-1.53944im 5.59917-4.04078im
 0.658702+6.55594im 1.23211+10.0951im 3.90816+2.07009im -2.79509+2.73057im 6.50781-2.18167im -2.71237-4.55771im
 -0.64458-4.56344im 4.39054-8.0552im -4.49956-3.66493im -1.37505-2.99915im -6.16341-6.68374im -2.62665-4.84961im
  5.83408+0.0im -1.513+5.81849im -9.69012-0.0557974im -6.16479+0.0im -9.69012+0.0557974im -1.513-5.81849im

```

(1,1), (1,4), and (4,1) are purely real. Furthermore (2, : ) is has the hermitian symmetry.

However, if you would naively cut the highest frequency of a 7x7 array to get a 6x6 array, you don’t have the symmetry:

```julia
julia> rfft(randn((7,7)))
4×7 Array{Complex{Float64},2}:
  3.66183+0.0im -1.02749+1.94808im 1.51224-14.965im -1.15603+6.66689im -1.15603-6.66689im 1.51224+14.965im -1.02749-1.94808im
 -1.03278+4.97944im -0.93996+8.06737im -6.88103+3.68525im 8.17271-4.17101im 4.03284+5.643im -0.919963-2.10859im -6.10537-1.09231im
 -1.03219+5.51279im 6.46411-11.5258im -0.71596-6.65423im -1.37923-10.1864im 0.996462+1.32339im -2.89172-1.68502im 4.21725-5.09607im
  0.34954+4.8042im 4.59035-4.21084im 0.637454-1.86879im 11.3116+6.97052im 5.91684+7.88591im 7.2835-2.23067im 1.83935+3.89034im

```

You are not providing the right symmetry `rfft` expects.

Let’s look at this, to downsample an odd sized 2D array to an even sized

```julia
function ds_2D(arr)
           a = rfft(arr)
           a = a[:, 1:end-1]
           a[:, end] .*= 0.5
           a[end, :] .*= 0.5
           a[end, end] *= 2
       
           return irfft(a, size(arr)[1]-1)
       end

ulia> rfft(x)
4×7 Array{Complex{Float64},2}:
  5.19375+0.0im -1.22145-3.05277im 1.99423+3.27891im -5.63567+2.6869im -5.63567-2.6869im 1.99423-3.27891im -1.22145+3.05277im
 0.468664+2.7739im -4.30717-2.67216im -1.11728+2.58736im -2.3358+3.52392im 0.419257-5.8995im -12.4438-5.26935im 0.413576+2.83563im
 -5.94199+4.45657im 4.97952+2.69954im -5.49123+7.09337im -1.05639+5.62915im -11.0779-6.8679im 7.01256+5.8489im 4.35303-3.09863im
 -4.61229-6.20257im -0.437476-8.4717im -3.10097-2.61191im 3.93253+5.40222im 1.05266-3.04878im 0.00995383-6.20056im 3.80061-8.75746im

julia> ds_2D(x) |> rfft
4×6 Array{Complex{Float64},2}:
  5.19375+0.0im -0.112167-0.706655im -1.82072+2.98291im -5.63567+0.0im -1.82072-2.98291im -0.112167+0.706655im
 0.468664+2.7739im -4.30717-2.67216im -1.11728+2.58736im -2.3358+3.52392im 0.419257-5.8995im -6.2219-2.63467im
 -5.94199+4.45657im 4.97952+2.69954im -5.49123+7.09337im -1.05639+5.62915im -11.0779-6.8679im 3.50628+2.92445im
 -2.30615+0.0im -0.106881-0.567786im -0.512076+0.109217im 1.96626+0.0im -0.512076-0.109217im -0.106881+0.567786im

```

You would think that the `ds_2D(x) |> rfft` is just a subset (maybe the highest frequency factor 0.5) of `rfft(x)` because we cut it in that way. But the first frequency line is completely different! And this will affect the 2D array because the artifacts are distributed over a whole frequency range in the first dimension.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 9, 2021, 12:08pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/53 "2021-01-09T12:08:00Z")

</div>

Why do you insist on `rfft()`? Just do regular `fft()`.  
As I wrote, for real signal add to my code `real()` to the Nyquist Sample.

Another way is do downsample by upsample.  
If your upsample works and validated, just do the trick I posted in [Sinc Interpolation based on FFT - #27 by RoyiAvital](https://discourse.julialang.org/t//52512/27).

Pay attention that while no LPF was applied this is the correct interpolation in the meaning those are the values of the “continuous” signal at those points.  
To apply Downsamplign in the meaning of estimating the value as the signal was sampled at lower rate you need to follow what I did on the DFT.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2021, 2:34pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/54 "2021-01-09T14:34:46Z")

</div>

> [@antoine-levitt](#):
>
> Not in this very particular case (downsampling from odd to even size). Of course the resulting imaginary part should be small (or your downsampling will be bad anyway and you should have filtered it before) but not of order of the machine precision.

I would still say you are doing it wrong if you are doing FFT-based downsampling and your algorithm doesn’t produce exactly conjugate-symmetric data in the frequency domain. Probably you have unbalanced handling of the Nyquist element(s).

By taking the real part of the output, you are throwing away the error you made, but as long as this error was in a small number of frequency components then it doesn’t affect the answer by much.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [January 9, 2021, 7:17pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/55 "2021-01-09T19:17:54Z")

</div>

My point is that taking the real part of the output is exactly equivalent to a proper handling of the nyquist elements, while being much simpler to write Anda lot less error prone. (Of course the proper way of doing things is exploiting the real FFTs, but I can never be bothered.)

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 9, 2021, 10:17pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/56 "2021-01-09T22:17:26Z")

</div>

> [@RoyiAvital](#):
>
> Pay attention that while no LPF was applied this is the correct interpolation in the meaning those are the values of the “continuous” signal at those points.  
> To apply Downsamplign in the meaning of estimating the value as the signal was sampled at lower rate you need to follow what I did on the DFT.

Yeah I agree on that.

> [@antoine-levitt](#):
>
> My point is that taking the real part of the output is exactly equivalent to a proper handling of the nyquist elements

For upsampling it’s reasonable, I agree on that.  
But for downsampling: I still don’t understand what’s the proper way? In 1D it is clear (multiplication of the Nyquist sample by 0.5). But for 2 or N dimensions, it isn’t that simple as I pointed out above, the symmetry handling is tricky.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 9, 2021, 10:34pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/57 "2021-01-09T22:34:12Z")

</div>

Am I wrong or did you manage to do the Upsampling for _Multi Dimensional Arrays_?  
If so, the logic of handling the Nyquist Sample for downsampling is exactly the same.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 9, 2021, 10:38pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/58 "2021-01-09T22:38:24Z")

</div>

I did, yes.

How can it be the same?  
For upsampling you append slices with hermitian symmetry to the other boundary, right?

For downsampling you naively would take the frequencies you want to keep, but then you need to fix the symmetry in the “Nyquist slices”?

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 9, 2021, 10:48pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/59 "2021-01-09T22:48:07Z")

</div>

The rationale of the adjustment for upsampling / downsampling is the same. The only difference is the indices.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [January 9, 2021, 10:52pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/60 "2021-01-09T22:52:19Z")

</div>

Sorry, I don’t really get what you mean by this 😕

[Previous page](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512.md?page=2)

[Next page](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512.md?page=4)