Sinc Interpolation based on FFT

No, at the Nyquist frequency it simply adds the aliased amplitudes like at every other frequency. There’s no “averaging”.

For example, for the IFFT of a “delta”-function frequency content, i.e. a pure sinusoid at a frequency of 2\pi*1/n (not the Nyquist frequency) with an amplitude of 1, it outputs an amplitude of 1. If we add an alias at a frequency of 2\pi*(1+n)/n then it outputs an amplitude of 2:

julia> using FFTW

julia> n = 8
8

julia> ifft( cis.(-2π * 1 * (0:n-1) / n) )
8-element Vector{ComplexF64}:
 -4.3063660604970826e-17 - 1.4302971815274583e-18im
                     1.0 + 1.244288779022145e-16im
   2.918587279715637e-17 - 1.5308084989341915e-17im
                     0.0 - 1.5308084989341915e-17im
  1.2447490626287002e-17 - 2.9185872797156375e-17im
 -1.1102230246251565e-16 - 3.2580367966163016e-17im
  1.4302971815274553e-18 - 1.5308084989341915e-17im
                     0.0 - 1.5308084989341918e-17im

julia> ifft( cis.(-2π * 1 * (0:n-1) / n) + cis.(-2π * (1+n) * (0:n-1) / n))
8-element Vector{ComplexF64}:
 -2.0859200112467696e-16 - 4.0288103043407935e-16im
                     2.0 + 9.560533348749494e-16im
  2.3634757674030587e-16 - 3.7512545481845044e-16im
                     0.0 - 3.3587314335135604e-16im
  -5.416589085122239e-16 + 9.671933064724106e-17im
                     0.0 + 4.065209743356281e-16im
    5.13903332896595e-16 + 6.896375503161215e-17im
                     0.0 - 4.1437776628554474e-16im

And exactly the same thing happens at the Nyquist frequency \pm \pi:

julia> ifft( cis.(-2π * (n/2) * (0:n-1) / n) )
8-element Vector{ComplexF64}:
                     0.0 - 6.123233995736767e-17im
 -2.5363265666181693e-17 - 6.123233995736766e-17im
  -6.123233995736764e-17 - 6.123233995736766e-17im
 -1.4782794558091698e-16 - 6.123233995736766e-17im
                     1.0 + 4.286263797015736e-16im
  1.4782794558091698e-16 - 6.123233995736766e-17im
   6.123233995736764e-17 - 6.123233995736766e-17im
  2.5363265666181693e-17 - 6.123233995736766e-17im

julia> ifft( cis.(-2π * (n/2) * (0:n-1) / n) + cis.(+2π * (n/2) * (0:n-1) / n) )

8-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 2.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

julia> ifft( 2cos.(-2π * (n/2) * (0:n-1) / n) )  # 2cos = same thing as adding alias
8-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 2.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

This shouldn’t be surprising since (a) the DFT definition does not have a different scale factor for the Nyquist term and (b) the aliased signal is identical at the sample points (this is the definition of “aliasing”) so it is literally impossible for the transform to detect how many aliased terms were summed in order to “average” them.

Maybe you are confused by the fact that if you put in a cosine (not a complex exponential) at a lower frequency, e.g. 2\pi/n, then a DFT is able to distinguish both the positive and negative frequency components as separate amplitudes:

julia> ifft( 2cos.(-2π * (1) * (0:n-1) / n) )
8-element Vector{ComplexF64}:
  -8.612732120994165e-17 + 0.0im
                     1.0 + 1.3973696289155642e-16im
  3.0616169978683824e-17 + 0.0im
 -1.1102230246251565e-16 + 1.7272282976821098e-17im
  2.4894981252574003e-17 + 0.0im
 -1.1102230246251565e-16 - 1.7272282976821098e-17im
  3.0616169978683824e-17 + 0.0im
                     1.0 - 1.3973696289155642e-16im

so the cosine seems to be “normalized differently” at the Nyquist frequency. This is just a consequence of the facts that the DFT basis functions are complex exponentials and a pure cosine consists of two complex exponentials, but at the Nyquist frequency those two complex exponentials are aliased.

2 Likes

You don’t need to do simulations, you can do this analytically. It simplifies a lot because (a) the DFT is linear and (b) the frequency components are orthogonal (when integrated over a common period), and the upshot is that you can analyze just the Nyquist term by itself:

In particular, take a function f(x)=a_{+}e^{i\pi x}+a_{-}e^{-i\pi x}, sampled at the Nyquist frequency f(n)=(a_{+}+a_{-})(-1)^{n} so we can only measure 2a=a_{+}+a_{-}. We reconstruct/interpolate it as a signal \tilde{f}(x)=2a\left[c_{+}e^{i\pi x}+c_{-}e^{-i\pi x}\right]. Question: what coefficients c_{\pm} minimize the expected mean-square error E[\int|f(x)-\tilde{f}(x)|^{2}dx] if a_{\pm} are i.i.d. random numbers with zero mean and some distribution (e.g. Gaussian)?

By explicit computation:

E\left[\frac{1}{2}\int_{0}^{2}|f(x)-\tilde{f}(x)|^{2}dx\right] \\ = \frac{1}{2} E\left[\int_{0}^{2}\left\{ |a_{+}-2ac_{+}|^{2}+|a_{-}-2ac_{-}|^{2}+\cancel{\#e^{2\pi ix}+\overline{\#}e^{-2\pi ix}}\right\} dx\right] \\ =E\left[|a_{+}-2ac_{+}|^{2}+|a_{-}-2ac_{-}|^{2}\right] \\ =E\left[|(1-c_{+})a_{+}|^{2}+|c_{+}a_{-}|^{2}+|(1-c_{-})a_{-}|^{2}+|c_{-}a_{+}|^{2}\right] \\ =E[|a_{\pm}|^{2}]\left(|1-c_{+}|^{2}+|c_{+}|^{2}+|1-c_{-}|^{2}+|c_{-}|^{2}\right) \\ =E[|a_{\pm}|^{2}]\left(1+2\left|\frac{1}{2}-c_{+}\right|^{2}+2\left|\frac{1}{2}-c_{-}\right|^{2}\right)

(using the facts that E[a_{+}\overline{a_{-}}]=0 and E[|a_{+}|^{2}]=E[|a_{-}|^{2}]), which is minimized for \boxed{c_{\pm}=\frac{1}{2}}.

That is, the optimal interpolant is \boxed{\tilde{f}(x)=a\left[e^{i\pi x}+e^{-i\pi x}\right]}, which corresponds to splitting the Nyquist amplitude 2a equally between the positive- and negative-frequency terms. This also coincides with the minimal mean-square slope interpolant from above, and has the nice property that it interpolates real signals f (a_- = \overline{a_+} \implies a purely real) with real interpolants \tilde{f}. (The same analysis also works, with the same optimal interpolant \tilde{f}, if we average over purely real signals f with random a_+=\overline{a_-}, since in that case we still have E[a_{+}\overline{a_{-}}]=E[a_+^2]=0 if the real and imaginary parts of a_+ are i.i.d. with zero mean.)

No \sqrt{2}.

8 Likes

Thanks for your nice examples (you may want to edit the result of abs. in line 3, block one). Yes, I do agree that it is the sum as I was previously thinking (Sinc Interpolation based on FFT - #72 by RainerHeintzmann) and how we first implemented FourierTools.jl. Yet what threw us off was the DFT of a delta-function placed at coordinate zero (or one in Julia notation).

julia> n=8; a = Base.setindex(zeros(n),1,1); abs.(fft(a))
8-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

Due to the aliasing, I would have expected a 2.0 at the 5th entry, but it is only a 1.0 (independently of even or odd n). I am still not sure how to understand this effect, but it must have a very simple explanation. It looks a bit like we have the choice to either preserve deltas or cosine functions correctly, when we downsample? But both the delta and the cosine are real and symmetric (leading to real and symmetric DFTs), so aliasing should act the same way on both. So what is going on here?

1 Like

Even Case (n = 8)

Just to make it concrete, let’s think in degrees rather than radians for a minute, you get the following for n = 8.

julia> [0:4; -3:-1]/8*360
8-element Vector{Float64}:
    0.0
   45.0
   90.0
  135.0
  180.0
 -135.0
  -90.0
  -45.0

julia> a = [0; 1; zeros(6)]
8-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

julia> fft(a)
8-element Vector{ComplexF64}:
                 1.0 + 0.0im
  0.7071067811865476 - 0.7071067811865476im
                 0.0 - 1.0im
 -0.7071067811865476 - 0.7071067811865476im
                -1.0 + 0.0im # 5th component
 -0.7071067811865476 + 0.7071067811865476im
                 0.0 + 1.0im
  0.7071067811865476 + 0.7071067811865476im

julia> fft(rand(8))
8-element Vector{ComplexF64}:
   3.1634222434593813 + 0.0im
   0.4951192018697982 + 0.4588330007918265im
 -0.11377277904888605 + 0.6583452797979255im
   0.4220473665291581 + 0.7562936096399402im
   -0.693998155489765 + 0.0im   # 5th component
   0.4220473665291581 - 0.7562936096399402im
 -0.11377277904888605 - 0.6583452797979255im
   0.4951192018697982 - 0.4588330007918265im

I think your confusion is that there should be a +180 and -180 degrees and thus the amplitude of that coefficient should be doubled. This is not the case because +180 and -180 are the same. If we shift the delta distribution over by one discrete unit, we see the 5th coefficient is real. In fact that will always be the case for a real input.

Odd Case (n = 9)

In the odd case +160 and -160 degrees are not the same and thus have distinct coefficients that are complex conjugates of each other.

julia> [0:4; -4:-1]/9*360
9-element Vector{Float64}:
    0.0
   40.0
   80.0
  120.0
  160.0
 -160.0
 -120.0
  -80.0
  -40.0

julia> a = zeros(9); a[2] = 1; fft(a)
9-element Vector{ComplexF64}:
                 1.0 + 0.0im
   0.766044443118978 - 0.6427876096865394im
 0.17364817766693036 - 0.984807753012208im
                -0.5 - 0.8660254037844386im
 -0.9396926207859084 - 0.3420201433256687im # 5th component (+160 degrees)
 -0.9396926207859084 + 0.3420201433256687im # 6th component (-160 degrees)
                -0.5 + 0.8660254037844386im
 0.17364817766693036 + 0.984807753012208im
   0.766044443118978 + 0.6427876096865394im

julia> fft(rand(9))
9-element Vector{ComplexF64}:
    5.806253282273899 + 0.0im
 0.027575433736157784 + 0.46275887781189695im
  -0.2594352504260514 + 0.06314184340863886im
   0.7252346656329784 - 0.8597235619145158im
  -1.3978344355575814 - 0.3593950913217314im # 5th component
  -1.3978344355575814 + 0.3593950913217314im # 6th component
   0.7252346656329784 + 0.8597235619145158im
  -0.2594352504260514 - 0.06314184340863886im
 0.027575433736157784 - 0.46275887781189695im

Video

I recommend picking up Grant Sanderson’s video at minute 13 where he gives a nice visualization of what is happening on the unit circle.

You may also be interested in his more casual video on the topic.

1 Like

Discrete/sampled signals are assumed to be periodic and band limited by the FFT technique.

An impulse at time zero is assumed to be repeating every N*dt samples and its FFT, from the formula definition, will be constant for all frequencies (in agreement with theory for the continuous case).

On the other hand, if a cosine signal is truncated in the time domain such that it does not have continuous integral periods every N*dt samples, then there will be additional frequency components (spectral leakage effect). In other words, the truncation of a periodic function over an interval that is not multiple of the period results in a sharp discontinuity in the time domain and in side-lobes in the frequency domain.

NB (a reminder for those who need it, including me):
In the continuous case, the Fourier Transform of c(t) = Acos(2pi f0 t) has two impulses with amplitude A/2 at f=+/-f0. If we let f0 to be 0, then c(t) = A (constant) and in the Fourier domain we get one single impulse of amplitude A at f = 0.

3 Likes

Discrete/sampled signals are assumed to be periodic and band limited by the FFT technique.

An impulse at time zero is assumed to be repeating every N*dt samples and its FFT, from the formula definition, will be constant for all frequencies (in agreement with theory for the continuous case).

On the other hand, if a cosine signal is truncated in the time domain, i.e., if it is not continuous every N*dt samples at the start/end of the window of analysis, then there will be additional frequency components (spectral leakage effect).

Yes, indeed. However here both examples are nicely repeated (the delta comb as well as the cosine) without any jumps.
Both examples are undersampled by the even-sized FFT. Both are real at the Nyquist frequency.
Yet one undersampling can be reverted easily, the other not. The problem is that the delta has lots of more undersampled frequencies whereas the cosine does not. Recovering the delta lives from the fake impression that we can recover also the higher frequencies but this is of course not true.
If you first band-limit the delta, than that band-limited signal behaves as you would expect.

The example and statement that follows are extremely puzzling:

Why one time series built with with one impulse at t=0 and a total of 8 samples would be aliased? And if it was why would it be aliased only in position 5 of the FFT and not in additional positions?
Probably this is being taken out of its proper context and therefore would appreciate if you could elaborate on the rationale here. Thank you.

Indeed it so puzzling which is why I posted it.
The cosine example posted before shown the factor of 2.0 for the aliased frequency at Nyquist sampling due to aliasing. Naive thinking made me expect the same for the delta peak also containing such a frequency which is aliased.
Yet it is not seen.
Remember that this “one delta pulse” really is a series of periodic delta pulses that do contain the frequency that is aliased.

Why is zero padding and Fourier transform the same as sinc interpolation?

Padding with M zeros, the interpolated frequency component k'\in [0, N+M-1] is

Y_{k'} = \frac{1}{N}\sum_{n=0}^{N+M-1} y_ne^{-2\pi i\frac{k'}{N+M}n} = \frac{1}{N}\sum_{n=0}^{N-1} y_ne^{-2\pi i\frac{k'}{N+M}n}

Plugging the FT y_n = \sum_{k=0}^{N-1} Y_k e^{2\pi i \frac{kn}{N}} into the previous expression, I get

Y_{k'} = \frac{1}{N}\sum_{k=0}^{N-1} Y_k \sum_{n=0}^{N-1} e^{2\pi i\left(\frac{k}{N} - \frac{k'}{N+M}\right)n}

The summation over n gives (using \sum_{n=0}^{N-1} q^n = \frac{1-q^N}{1-q})

K(k,k')\equiv \frac{1- \exp\left(2\pi i\left(\frac{k}{N} - \frac{k'}{N+M}\right)N\right)}{1 - \exp\left(2\pi i\left(\frac{k}{N} - \frac{k'}{N+M}\right)\right)} = \frac{\sin\left(\pi \left(\frac{k}{N} - \frac{k'}{N+M}\right)N\right)}{\sin\left(\pi \left(\frac{k}{N} - \frac{k'}{N+M}\right)\right)}e^{\pi i\left(\frac{k}{N} - \frac{k'}{N+M}\right)(N-1)}

and

Y_{k'} = \frac{1}{N}\sum_{k=0}^{N-1} Y_k\cdot K(k,k')

Where is the sinc function?

Note that in this case we mean the periodic sinc function also referred to as the Dirichlet function or kernel.

https://math.stackexchange.com/questions/2306017/sinc-function-vs-dirichlet-kernel/2306055

1 Like

Thank you for the pointer, but I still do not get it how this is related to the discrete Fourier transform. How is the Dirichlet function D_n(x) used in the discrete Fourier transform?

You could probably post the simple case where M = N. Also, over the range of indices considered, k' = k, and so the formulas will simplify greatly and it will be easier to reason.

It comes from the DFT of a rectangular window. Zero-padding can be thought of as multiplying a longer signal (which you don’t have) by a square window (resulting in the data you have), and hence by the convolution theorem this takes the spectrum of the longer signal and convolves it with a Dirichlet kernel, which is a form of trigonometric interpolation.

4 Likes