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

<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 10, 2021, 8:40am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/61 "2021-01-10T08:40:24Z")

</div>

> [@roflmaostc](#):
>
> 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.

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?

---

<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 10, 2021, 10:39am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/62 "2021-01-10T10:39:58Z")

</div>

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

```julia
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

```

```julia
# 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`:

```julia
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](https://dsp.stackexchange.com/a/72434).

---

<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 10, 2021, 2:49pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/63 "2021-01-10T14:49:57Z")

</div>

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.

---

<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 10, 2021, 3:10pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/64 "2021-01-10T15:10:04Z")

</div>

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](https://github.com/roflmaostc/FFTInterpolations.jl)  
And I leave it like this for now.

Thanks for your input!

---

<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 10, 2021, 7:29pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/65 "2021-01-10T19:29:54Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [January 16, 2021, 7:52pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/66 "2021-01-16T19:52:57Z")

</div>

@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](https://github.com/roflmaostc/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:

 ![1d-sinc_interpolation](https://global.discourse-cdn.com/julialang/original/3X/a/9/a95ce240dd70178fc143ddcd5ac84cad92bf6cdc.png)

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](https://github.com/roflmaostc/FFTResampling.jl):

 ![2d-sinc_interpolation](https://global.discourse-cdn.com/julialang/original/3X/5/0/50eed6b66bf9860910d2e95b165e15cd10dbe449.png)

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.

---

<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 17, 2021, 1:00pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/67 "2021-01-17T13:00:27Z")

</div>

I was also considering to rename it and you gave men the final push 😉  
[FFTResampling.jl](https://github.com/roflmaostc/FFTResampling.jl)

---

<div class="post-metadata">

**Author:** ![Arrigo\_Benedetti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arrigo_benedetti/32/25545_2.png) [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)\
**Post date:** [February 8, 2021, 4:47pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/68 "2021-02-08T16:47:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [February 22, 2021, 12:36pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/69 "2021-02-22T12:36:26Z")

</div>

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

```julia
resample = N_new > N_old ? sinc_interpolate : downsample
A_new = resample(A_old, N_new) 

```

---

<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:** [February 22, 2021, 12:43pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/70 "2021-02-22T12:43:19Z")

</div>

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

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

---

<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:** [February 23, 2021, 7:43pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/71 "2021-02-23T19:43:13Z")

</div>

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

```julia
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](https://github.com/roflmaostc/FFTResampling.jl/blob/03a6253f68cfc4a717bd006ec3a2e28fcb0e5159/src/sinc_interpolation.jl#L6)

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.

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [March 30, 2021, 7:05pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/72 "2021-03-30T19:05:32Z")

</div>

We now implemented FourierTools.jl a first version of which will appear soon. I want to use this discussion to outline how a proper padding and cutting in Fourier space (also mixed versions) needs to be done.

1. We want: Reversibility for cases where we do not reduce sizes.
2. We need to not change the integrated power spectrum

The latter part is crucial when it comes to dealing with even sized dimensions in FFT and RFFT of real valued data alike. The lowest (biggest negative) frequency stores in these cases the sum of amplitude of the negative highest frequency PLUS the corresponding value at the aliased positive highest frequency (e.g. [-kxmax, ky, kz] + [+kxmax, ky, kz]). Note that this value is generally NOT real except for special choices of ky and kz.

When changing the region of interest in FFT- (or RFFT-) space the algorithm consists of two steps:

1. `fix_dim_before(Array, dim)`
2. `fix_dim_after(Array, dim)`  
Both of these routines need to be run for each dimension consecutively. The first before the size change and the second after it.

`fix_dim_before(Array, dim)` is only active for a dimension if it is even-sized and enlarged  
In this case the highest frequency subslice(dim) needs to be halfed in value and this half also copied (without conjugation) to the equivalent positive frequency (e.g. [+kxmax, ky, kz]).

`fix_dim_after(Array, dim)` is only active for a dimension if it is even-sized and reduced  
Here we need to copy the lowest frequency which is just about dropped (e.g. [+kxmax, ky, kz]) and sum it to the other side (here [-kxmax, ky, kz]))

We implemented this for array views (overloading `getindex`, such that no extra memory allocation is needed.  
A remark on just taking the real part and ignoring such (rather complicated) actions:  
An even lowest frequency in FFTs (and RFFTs) store both, the positive and negative highest frequencies. Going back to real-space after padding in Fourier-space without copying will then drop some of the energy at the moment when the real part is taken. When reverting this situation (cutting away the previously introduced zeros) it is clear to see that this real-space energy is still lost.

An interesting consequence is that also some other functions need modifications on which we are currently working: For example: multiplication with a complex exponential needs a modification to the lowest frequency along an even dim.  
It seems to be that also calculating the Energy or Power spectrum needs a modification to such frequencies: sum((a+a)^2) !== 2\*sum(a^2)  
The should be fairly easy to check via a comparison to the real-space integral of the square of the signal.  
If array sizes are fairly large, all of these effects are small, but wouldn’t it be nice to get this correct?  
This is what `FourierTools.jl` will be trying to achieve.

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [April 1, 2021, 10:14am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/73 "2021-04-01T10:14:30Z")

</div>

I have to correct myself: The lowest frequency in an FFT does NOT represent the sum of this frequency content and the corresponding (aliased) frequency from the positive max frequency, but it corresponds to the AVERAGE. This is easily seem by performing an FFT of a delta-shaped distribution at zero coordinate: `using FFTW, IndexFunArrays; fft(delta((4,4),offset=(1,1)))`.  
This does make a significant difference for the necessary operations and lends validity to the real-casting in real space. Yet, I think this only works in 1D.  
Discussing with @roflmaostc we concluded that the essential problem with even-sized arrays is the following ambiguity: The highest frequency may represent a cosine-shaped wave with phase 0 or any phase shifted wave with a much higher amplitude but a different phase. As opposed to odd-sized arrays: Not only can the FFT highest frequency not distinguish between these cases but also the real-space array itself is ambiguous in this respect. In Fourier space, the lowest (highest negative) frequency is really the sum of two (unknown) phasors with unknown relative phase angle and corresponding magnitude such that the average yields the stated real part.  
This means that there is really no unambiguous “correct” way to disentangle this, but one can argue about different ways to do this.

1. Assuming the (fundamentally unknown) components to have real-only phase.  
The consequence is that on average (randomly) distributed phases, we are overestimating the power of this frequency. Yet the advantage is that centrosymmetric distributions (known to be real in Fourier space) are fine. Splitting is (opposed to what I wrote before) not required to halve the value, but just to copy it upon padding in Fourier space.
2. We keep real-valuedness but mulitply this frequencies by sqrt(2) upon splitting. This assumes that these frequency did realy have an unkown randomized phase before they were averaged together.  
Since we have no chance to correct the (unkown) phase, we can at least try to correct the statistical effect on the magnitude the averaging was having. Advantage: The power density after padding (or even before padding) has ON AVERAGE a better representation. Disadvantage: For centrosymmetric structures we get a too strong powerspectrum. Disadvantage: Upon clipping the array (back) what do we do then? To be reversible in the operation it sounds like we would need to devide be sqrt(2) after averaging the corresponding pixels. But here we actually do know the phases and simple averaging is the best result.
3. We set the highest frequency to zero before padding. This is a loss of data which probably should be avoided if possible.

Taken together it looks like option one is the way to go.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [April 2, 2021, 2:40pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/74 "2021-04-02T14:40:16Z")

</div>

@stevengj has expertise on this topic.

---

<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:** [April 2, 2021, 2:55pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/75 "2021-04-02T14:55:04Z")

</div>

> [@RainerHeintzmann](#):
>
> In Fourier space, the lowest (highest negative) frequency is really the sum of two (unknown) phasors with unknown relative phase angle and corresponding magnitude such that the average yields the stated real part.

Yes, the Nyquist component is ambiguous in this sense, but there are ways to disambiguate. For example, for a band-limited signal there is a unique choice that minimizes the mean-square slope of the interpolated function, and the same choice also has the nice property of giving a real interpolant from real inputs. I have some notes on this: [https://math.mit.edu/~stevenj/fft-deriv.pdf](https://math.mit.edu/~stevenj/fft-deriv.pdf)

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [April 2, 2021, 4:40pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/76 "2021-04-02T16:40:09Z")

</div>

I am not sure, I understand. I think this ambiguity is pretty fundamental. Let me give an example:  
A cosine curve at the limit with even array size yields maximum contrast, but a sine curve is entirely invisible as it yields only zeros. If we now sample a function with a phase in between the two extreme cases, we would always interprete this as a cosine-only. Maybe it is the best we can do, but it is unfortunately not correct. Of course this is NOT a problem of the FFT but a problem of just about incorrect sampling as we are AT the Nyquist limit rather than sampling above it. @stevengj do you agree?

---

<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:** [April 2, 2021, 5:46pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/77 "2021-04-02T17:46:57Z")

</div>

> [@RainerHeintzmann](#):
>
> @stevengj do you agree?

This is just an example of [aliasing](https://en.wikipedia.org/wiki/Aliasing). Yes, **looking only at the samples** n it is impossible to distinguish between 2a\cos(\pi n) = ae^{i\pi n} + ae^{-i\pi n} and 2a\cos(\pi n) + 2b\sin(\pi n) = (a-ib)e^{i\pi n} + (a+ib)e^{-i\pi n} for any b, because the sine term vanishes on all the samples (equivalently, e^{\pm i \pi n} are aliased). You _must_ assume something about what happens in between the samples to determine b. For example, if you want the interpolant that minimizes the mean-square slope you get b=0, since the mean-square slope is:

\frac{1}{2} \int\_{-1}^{+1} \left| -2\pi a\sin(\pi x) + 2\pi b\cos(\pi x) \right|^2 dx \\ = 2\pi^2 \int\_{-1}^{+1} \left( |a|^2 \sin^2(\pi x) + |b|^2 \cos^2(\pi x) - \Re[a^\* b] \sin(2\pi x) \right) dx \\ = \pi^2 \left( |a|^2 + |b|^2 \right)

which is obviously minimized for b=0, whereas a is fixed by the sample observations.

Interpolating inherently involves choices, because there are infinitely many functions that interpolate between a discrete set of samples. Trigonometric interpolation is no different. Fortunately, there are reasonable criteria (like bandlimited interpolants with minimal mean-square slope) that lead to a simple choice of interpolant.

> [@RainerHeintzmann](#):
>
> Of course this is NOT a problem of the FFT but a problem of just about incorrect sampling as we are AT the Nyquist limit rather than sampling above it.

In a sense, yes — if you sample a DFT exactly at the Nyquist limit then the bandlimited assumption is not sufficient to determine the interpolant. You need an additional criterion like the minimal-slope condition above. But I don’t see this as a big problem; after all, the band-limited assumption itself is usually just an approximation.

(For the classic form of the Nyquist–Shannon sampling theorem, you have infinitely many samples — a [DTFT](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform), not a DFT — giving a continuous set of frequencies in [-\pi,\pi] and corresponding amplitudes. In this case, the Nyquist frequency \pm \pi is a set of measure zero and doesn’t contribute anything to the interpolated signal as long as the Fourier transform is [regular](https://www.sciencedirect.com/topics/mathematics/regular-distribution) (no delta functions). Hence the bandlimited assumption is sufficient by itself to determine a unique interpolant/reconstruction.)

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [April 2, 2021, 8:32pm UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/78 "2021-04-02T20:32:19Z")

</div>

Thanks for the detailed explanation. As I am in the field of optics, we are in the lucky position that for the ideal (i.e. noise-free) signal the band limit is indeed fulfilled.  
On a more practical side though: When it comes to implementing the Fourier-padding as we are attempting in our FourierTools package right now, we did implement it under the assumption of the highest (1D) frequency being real, which corresponds to the case you mention (minimal curvature) as far as I can see. Yet, if I understand it correctly this frequency in a DFFT is NOT just the aliasing, which would be a sum, but it is an average. But if we have very informative signals, they should be close to random and this would probably mean a more or less random phase. And this would mean that “on average” we would be slightly better off to assume this value to be representing a by sqrt(2) reduced version. Agreed?

I am not insisting to implement it like this as many practical signal are symmetric (i.e. real for this frequency). On the other hand, if we have a symmetric signal that is randomly offset by subpixel positions: There you go. This frequency may be on average better represented by the sqrt(2) assumption.

---

<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:** [April 3, 2021, 2:16am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/79 "2021-04-03T02:16:18Z")

</div>

> [@RainerHeintzmann](#):
>
> Yet, if I understand it correctly this frequency in a DFFT is NOT just the aliasing, which would be a sum, but it is an average.

Aliasing means you can’t distinguish e^{\pm i\pi x}, i.e. frequencies \pm \pi, and can only obtain the sum of the two coefficients, which is exactly what is going on here as I explained above (you get 2a = (a-ib) + (a+ib)).

> [@RainerHeintzmann](#):
>
> But if we have very informative signals, they should be close to random and this would probably mean a more or less random phase. And this would mean that “on average” we would be slightly better off to assume this value to be representing a by sqrt(2) reduced version. Agreed?

I don’t know what this means. What is an “informative signal”? Are your signals actually periodic, or are you using the DFT as a windowed approximation for the DTFT? Is this measured data that is sampling a physical continuous-time signal that you are trying to reconstruct, or…

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [April 5, 2021, 11:28am UTC](https://discourse.julialang.org/t/sinc-interpolation-based-on-fft/52512/80 "2021-04-05T11:28:05Z")

</div>

> > ![](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/40/19726_2.png)RainerHeintzmann:
> > 
> > Yet, if I understand it correctly this frequency in a DFFT is NOT just the aliasing, which would be a sum, but it is an average.
> 
> Aliasing means you can’t distinguish e^{\pm i\pi x}, i.e. frequencies \pm \pi, and can only obtain the sum of the two coefficients, which is exactly what is going on here as I explained above (you get 2a = (a-ib) + (a+ib)).

I know what aliasing means. The point I am making is the instead 2|a|, which is what I would have expected at first, the DFFT as implemented yield only |a|. You can check for yourself by Fouriertransforming a delta function on an even sized grid.

> > ![](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/40/19726_2.png)RainerHeintzmann:
> > 
> > But if we have very informative signals, they should be close to random and this would probably mean a more or less random phase. And this would mean that “on average” we would be slightly better off to assume this value to be representing a by sqrt(2) reduced version. Agreed?
> 
> I don’t know what this means. What is an “informative signal”? Are your signals actually periodic, or are you using the DFT as a windowed approximation for the DTFT? Is this measured data that is sampling a physical continuous-time signal that you are trying to reconstruct, or…

I agree, “informative” is not really an established scientific term. What I meant is that the signal looks close do noise (== cannot be compressed == informative) as opposed to something that is for example known to be real and symmetric. Let’s assume that the signal is band-limited including only the highest frequency in question but nothing beyond. Maybe I can do a simulation of 100s of noisy signals which are filtered, to demonstrate how you probably get get a better estimate by the sqrt(2) scaling at the Nyquist sampling limit than by the mean (or sum as you mentioned).

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

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