Is the discrete Gaussian kernel an eigenfunction of the DFT?

So the Gaussian function is an eigenfunction of the Fourier transform because it transforms into itself, right?

But this isn't true for the sampled Gaussian in the DFT because the tails of the function are truncated, right?

Wikipedia describes a discrete Gaussian kernel here and here (solid lines), which is different from the discretely-sampled Gaussian (dashed lines):

the discrete counterpart of the continuous Gaussian in that it is the solution to the discrete diffusion equation (discrete space, continuous time), just as the continuous Gaussian is the solution to the continuous diffusion equation Does that mean it also DFT transforms exactly into itself? If not, is there a similar Gaussian-like function that does?

• Not necessarily (or ever?) Gaussian-like, but see @robertbristow-johnson's comments to Interpolation of magnitude of discrete Fourier transform (DFT). Aug 17 '20 at 6:43
• @OlliNiemitalo Olli, using a Gaussian Eignenvector makes the "frequency humps" in the spectrum the same shape no matter where the frequency is located relative to bins, thus interpolation (with a single tone signal) will yield an exact answer. That's what is exciting for me about finding the field where they grow. The other piece of news (working on proving this) is that there is only one minimal width eigenvector per N (per type) and the rest are composed from this. Aug 17 '20 at 12:40
• @OlliNiemitalo (avoiding the chat warning) I solved the two sample peak problem by multiplying the signal by a pure complex tone to slide the spectrum peak up from DC to matching the signals peak between the bins. The spectrum then has a two sample peak that matches. Reversing that means the eigenvector is going to need to be spun up some complex tone making it not real. Negative values are expected in the higher order eigenvectors. There's a simple demo program here dsprelated.com/thread/12089/has-any-one-seen-this-window-family (and the throwing of a gauntlet). Aug 19 '20 at 18:13
• @OlliNiemitalo It should also be possible to add the W_L member with shifted copy of itself to create a single sample peak, but this widens the stance a little bit. No problem. Likewise, it should be possible to add the vectors in some binomial weighted manner to create a wider eigenvector for the same N. I'm still trying to get this figured properly. That's the larger puzzle here. Aug 19 '20 at 18:40
• @OlliNiemitalo They are definitionally the same, so yes, those are pure eigenvectors of your rotated DFT with the right b values. The code confirms this. Aug 20 '20 at 13:09

Since the DFT is representable by multiplication with the Fourier matrix, your question is equivalent to asking what the eigenvectors of the Fourier matrix are.

Actually, Wikipedia provides the answer (http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Eigenvalues_and_eigenvectors).

However, since the eigenvalues ($1, -1, i, -i$) are not simple, the eigenvectors are not unique (i.e. linear combinations are also eigenvectors). Also no simple closed formula exists.

One formula for an eigenvector close to what you ask is provided by Wikipedia

$$F_m = \sum_{k = -\infty}^\infty \exp\left(-\frac{\pi\cdot(m+N\cdot k)^2}{N}\right) \quad\quad\quad m = 0,\ldots,N-1$$

So concluding, the Gaussian function itself is not an eigenvector, but an infinite sum of Gaussians. The infinite sum can probably be interpreted as equivalent to the discretization of the frequency and time domain when we going from the FT to the DFT. So it is not as simple, as just truncating the discrete Gaussian.

• Isnt an infinite sum of gaussians still a gaussian though? Aug 22 '13 at 15:43
• No, convolution of Gaussians is still a Gaussian. The sum is Gaussian only if they are the same position and width. This function here is actally one period of a discrete Gaussian pulse train. So it doesnt even look like a Gaussian. Aug 22 '13 at 15:55
• Ah, I see. In other words this sum is essentially a Gaussian train made up of Gaussians of same variance but different means? Aug 22 '13 at 15:59
• exactly. The means are spaced by exactly N, the length of the DFT. Aug 22 '13 at 16:01
• Ah, fascinating. One last thing, this is an infinite length vector, which means that the DFT matrix is also infinite in length, no? Aug 22 '13 at 16:03

I have made a tremendous amount of progress on this issue in the last few weeks.

I have discovered a previously unrecognized class of window functions. What is particularly amazing about these is that the Eigenvectors of the DFT seem to be regularly placed members in this family. This provides not only a numerical solution for calculating the exact minimal sized eigenvectors, but also a concise definition of the family which should be usable in the proof they are eigenvectors, but also in constructing new generators for the higher order discrete Hermite-Gaussian functions.

I am still working out the details of how all the larger family fits together and searching for difference equations to define the eigenvectors from initial conditions.

There will no sleep lost in this house if someone beats me to it.

This is the definition of the Zeroing Sine Windows:

$$W_L[n] = \prod_{l=1}^{L} \sin \left( \frac{n+l}{N}\pi \right)$$

This is the definition of Raised Sine Windows:

$$W_L[n] = \sin^L \left( \frac{n}{N}\pi \right)$$

The former is the "factorial" version, the latter the "power" version. They start out being similar.....

Update: The recursive function for generating a single member has been found:

\begin{aligned} W_L &= 0.0001 \\ W_L[n] &= \frac{\sin \left( \frac{n+L}{N}\pi \right)}{\sin \left( \frac{n}{N}\pi \right)} W_L[n-1] \\ \end{aligned}

It comes straight from the definition.

The constant doesn't matter as long as it isn't zero. The results will be proportionate.

With $$N$$ odd, and $$L = (N-1)/2$$, this will be the minimum width eigenvector of the DFT.

See my article for the even case fix.

For any doubters, here is a quick sample program. Only the odd N of the form $$4m+1$$ will have an eigenvector with a single sample at the peak. $$L=(N-1)/2$$ and the peak occurs at $$L/2$$. The other N's need a little massaging to put them into centered eigenvector form.

import numpy as np

#========================================================================
def main():

#---- Set parameters

m = 4

N = m * 4 + 1
L = ( N - 1 ) >> 1

#---- Build the Zeroing Sine Family member

w = np.zeros( N )

w = 2**(-L)

for n in range( 1 , N ):
w[n] = w[n-1] * np.sin( (n+L) * np.pi / N ) / np.sin( n * np.pi / N )
print( " %3d %12.8f" % ( n, w[n] ) )

#---- Shift the function to be zero centered

r = np.roll( w, N - (L >> 1) )

#---- Take the DFT and compare

sr = np.fft.fft( r ) / np.sqrt( N )

for n in range( N ):
if r[n] == 0 :
print( " %3d %12.8f %12.8f %12.8f     Undefined" % \
( n, r[n], sr[n].real, sr[n].imag ) )
else:
print( " %3d %12.8f %12.8f %12.8f  %12.8f" % \
( n, r[n], sr[n].real, sr[n].imag, sr[n].real / r[n] ) )

#========================================================================
main()

Here are the results.

1   0.02116787
2   0.05636062
3   0.09583753
4   0.11352307
5   0.09583753
6   0.05636062
7   0.02116787
8   0.00390625
9   0.00000000
10  -0.00000000
11   0.00000000
12  -0.00000000
13   0.00000000
14  -0.00000000
15   0.00000000
16  -0.00000000
0   0.11352307   0.11352307   0.00000000    1.00000000
1   0.09583753   0.09583753  -0.00000000    1.00000000
2   0.05636062   0.05636062  -0.00000000    1.00000000
3   0.02116787   0.02116787  -0.00000000    1.00000000
4   0.00390625   0.00390625   0.00000000    1.00000000
5   0.00000000   0.00000000   0.00000000   26.92919942
6  -0.00000000   0.00000000   0.00000000  -295.67062572
7   0.00000000   0.00000000   0.00000000  1209.76913813
8  -0.00000000   0.00000000   0.00000000  -826.54020309
9   0.00000000   0.00000000  -0.00000000  826.54020309
10  -0.00000000   0.00000000  -0.00000000  -1209.76913813
11   0.00000000   0.00000000  -0.00000000  295.67062572
12  -0.00000000   0.00000000  -0.00000000  -26.92919942
13   0.00390625   0.00390625  -0.00000000    1.00000000
14   0.02116787   0.02116787   0.00000000    1.00000000
15   0.05636062   0.05636062   0.00000000    1.00000000
16   0.09583753   0.09583753   0.00000000    1.00000000
• +1, they are the real thing Aug 18 '20 at 12:14
• What does the term discrete Hermite-Gaussian function refer to? Aug 19 '20 at 4:46
• @OlliNiemitalo They are "eigenfunctions of the DFT" (search term). They are the leading polynomials of the Gaussian to make the eigenvectors. My gut feeling is they correspond to using different $W_0$ functions in my formulation (like the "crazy uncle" even case) and they can be found by simple differences (with perhaps twists) of the lower level ones. Try clustering most the zeroes on the left side, but place the last one in the middle of the non-zero interval. I'm probably going to give it another day (maybe two) before I start experimenting with this to find out. Go Olli! Aug 19 '20 at 13:03
• Yes it is very annoying that they keep papers behind a paywall. Waste of useful information. You can get that paper for free from Researchgate. Many paper PDF files can be found by googling for the title + pdf. The centered DFT might be something like this one with $b=1/2$. If I were to dabble with the math, I'd do everything centered around time and frequency 0, to stay within real numbers. Aug 19 '20 at 13:35
• I could not find the google drive link in the 2015 paper, but I found it in Alexey Kuznetsov, Mateusz Kwaśnicki, 2017, Minimal Hermite-type eigenbasis of the discrete Fourier transform. Aug 22 '20 at 13:12

This answer complements @CedronDawg's answer which introduced this family of eigenvectors. More specifically, this answer presents three algorithms and a hybrid algorithm for generating for a given length $$N$$ the unique (except for normalization) non-negative eigenvector of discrete Fourier transform (DFT) with the maximal number of consecutive zeros. This is the normal DFT with no half-sample shifts in time or frequency.

Iterative algorithm, $$\sim O\big(\exp(N)\big)$$ time

Here is an iterative numerical algorithm to find the eigenvectors satisfying the said constraints. The algorithm principle is to repeat enforcing the defining constraints in each domain and averaging the results. In Octave:

function retval = constrain(x, N)
zerocount = N-floor((N+2)/4)*2-1;
if (zerocount > 0)
firstzero = floor((N+10)/4);
x(firstzero:firstzero+zerocount-1) = 0;
endif
x = abs(x);
retval = x / x(1);
endfunction

function retval = geteig_iterative(N, precision = 0, x = ones(1, N))
x = constrain(x, N);
good = x;
X = fft(x)/sqrt(N);
error = sumsq(x-X);
iter = 0;
while (max(abs(X-x)) > precision)
lasterror = error;
X = constrain(X, N);
temp = X;
x = (x + X)/2;
X = fft(x)/sqrt(N);
error = sumsq(x-X);
iter++;
if (error >= lasterror)
break;
endif
good = temp;
endwhile
printf("Iterations: %d\n", iter);
retval = good;
endfunction

N = 21; x = geteig_iterative(N); plot([0:N-1], x, "x", [0:N-1], real(fft(x)/sqrt(N)), "+");

The algorithm can be given a starting point vector. Otherwise it starts with a vector of ones. The parameter precision defines the desired precision as the maximum absolute difference between the calculated eigenvector and its DFT, for example precision = 0.000000000001, or the default precision = 0 to halt the algorithm when an iteration no longer improves the solution.

The variables firstzero and zerocount determine the location and the number of zeros. Their formulas are in the realm of empirical mathematics, but were also guided by experience from symbolically solving the problem for some small values of $$N$$, and partially by testing whether random-vector initialization changes the results or if zerocount can still be increased.

For $$N = 21$$, the result is obtained in 2605 iterations: Figure 1. The output of the iterative algorithm (✕) and its DFT (＋), for $$N = 21$$.

Unfortunately, the algorithm needs a number of iterations that is exponentially proportional to $$N$$, so it is not useful for large $$N$$, and it also fails to converge for some not that large values of $$N$$ like $$18$$ and $$22$$: Figure 2. Number of iterations needed for precision = 0.000000000001, as function of $$N$$.

Despite the algorithm's shortcomings at large $$N$$, the results at smaller $$N$$ are useful, showing clear structure suitable for algorithmic construction. The structure is best viewed by looking at the frequency-domain and other z-plane zeros of the non-zero portion of the unwrapped eigenvectors, in Octave by zplane(nonzeros(fftshift(geteig_iterative(N)))');: Figure 3. $$z$$-plane zeros (◯) of the non-negative subsequences of the unwrapped eigenvectors. Poles (✕) can be ignored in these Z-plane plots. Not too easy to see, but there is a double zero at $$z = -1$$ when $$\operatorname{mod}(N, 4) \in \{0, 2\}$$. $$N < 3$$ were not included, because it was ambiguous how to split the unwrapped eigenvector to obtain the non-zero subsequence.

The analysis shows a $$4$$-periodic pattern in the zero distributions, as function of $$N$$. Double zeros at $$z = -1$$ can be seen when $$\operatorname{mod}(N, 4) \in \{0, 2\}$$. The two first columns with $$\operatorname{mod}(N, 4) \in \{0, 1\}$$ show a regular spacing of the zeros on the unit circle, representing real zeros in frequency domain that correspond to the frequencies of the DFT bins. For $$\operatorname{mod}(N, 4) \in \{2, 3\}$$, there are two necessary but somewhat arbitrarily placed real zeros, one inside and one outside the unit circle, with their value uniquely defined for a given $$N$$, as confirmed by testing with random-vector initializations. Symmetry constraints require that the additional pair of zeros is a reciprocal pair, reducing the number of unknowns to just one, which is an easy problem to solve numerically, as will be done in the fast algorithms in the following sections.

The simplicity of the pattern of zeros seen with $$\operatorname{mod}(N, 4) = 1$$ is attractive. It may be attractive enough to choose a time and/or frequency-shifted variant of DFT depending on $$\operatorname{mod}(N, 4)$$, to ensure the same simplicity for all $$N$$. However, it may be important to include a zero-frequency bin, and as the question was about the usual kind of DFT, I will stick to that for this answer.

Fast convolutional algorithm, $$O\big(N\log^2(N)\big)$$ time

@CedronDawg's approach was to construct the eigenvectors as products of real-valued vectors that each created at most two zeros. Here is presented a fast algorithm that uses a similar principle, but instead of multiplying such vectors one by one, it takes a divide-and-conquer approach. In this convolutional algorithm, frequency-domain zeros are created by time-domain convolution. Convolution by conv is accelerated by fast Fourier transform (FFT), which gives a total time complexity of $$O\big(N\log^2(N)\big)$$.

The main work in the algorithm is to create a long train of frequency-domain zeros. Such work is split to creating and combining a pair of half-length trains, recursively, until the length of the train to create is only 3 or less. The time-domain sequence corresponding to such a short train can be expressed simply and directly. For odd-length trains the middle zero is created separately. Instead of creating half-trains individually, they are cloned and shifted to the appropriate place by complex sinusoid modulation. This way the recursions do not create a call tree but only a short call chain of length $$O\log(N)$$. In Octave:

function retval = getzeros(N, M)
if (M == 1)
retval = [1, -1];
elseif (M == 2)
retval = [1, -2*cos(pi/N), 1];
elseif (M == 3)
retval = [1, -1 - 2*cos(2*pi/N), 1 + 2*cos(2*pi/N), -1];
elseif (rem(M, 2) == 0)
half = getzeros(N, M/2);
modulator = exp(i*[((0:M/2)-M/4)*2*pi*M/4/N]);
retval = real(conv(conj(modulator).*half, modulator.*half));
else
half = getzeros(N, (M-1)/2);
modulator = exp(i*[((0:(M-1)/2)-(M-1)/4)*2*pi*(M+1)/4/N]);
retval = conv(real(conv(conj(modulator).*half, modulator.*half)), getzeros(N, 1));
endif
endfunction
#[h, w] = freqz(getzeros(10, 4), [], 500); plot(abs(h));

function retval = geteig_convolutional(N)
M = N-floor((N+2)/4)*2-1;
x = getzeros(N, M);
if (rem(M, 2) == 1)
x = conv(x, getzeros(N, 1));
endif
x .*= (-1).^[0:length(x)-1];
mid = (length(x) + 1)/2;
if (rem(N, 4) == 2 || rem(N, 4) == 3)
s = sum(x)/sqrt(N);
x = conv(x, [1, 2*(s - x(mid + 1))/(x(mid) - s), 1]);
endif
retval = circshift(horzcat(x/x(mid), zeros(1, N-length(x))), (rem(N,4) < 2)-mid);
endfunction

N = 512; x = geteig_convolutional(N); plot([0:N-1], x, "x", [0:N-1], real(fft(x)/sqrt(N)), "+");
plot([0:N-1], x - real(fft(x)/sqrt(N)), "x");

This plots: Figure 4. The output of the fast convolutional algorithm (✕) and its DFT (＋), for $$N = 512$$. Figure 5. The difference between the output of the fast convolutional algorithm and its DFT, for $$N = 512$$.

The convolutional algorithm begins to explode with numerical error for large values of $$N$$ roughly $$N > 750$$, I guess partially because of back-and-forth FFT – inverse FFT (IFFT). The multiplicative algorithm presented in the next section will be much more stable.

The code includes creating the reciprocal zeros needed for $$\operatorname{mod}(N, 4) \in \{2, 3\}$$ as seen in Fig. 3. For that, the value for the variable $$a$$ in convolution by $$[1, a, 1]$$, which creates the reciprocal zeros, was solved from an equation $$x[-1] + a\,x + x = X + a\,X + X$$. The equation requires that after convolution, the value at time domain origin (represented by the left side of the equation) equals the value at frequency domain origin (right side), with $$x$$ and $$X$$ representing the nascent eigenvector and its DFT before the convolution.

Fast multiplicative algorithm, $$O\big(N\log(N)\big)$$ time

Here is the multiplicative version of the fast algorithm. It creates zeros in the same domain in which it constructs the eigenvector. Intermediate results are stored in logarithmic scale to accommodate for the large dynamic range of the numbers. In Octave:

function retval = get_log2_zeros_multiplicative(N, M)
k = [0:N-1];
if (M == 1)
x = log2(abs(sin(pi*k/N)));
elseif (M == 2)
x = log2(abs(cos(pi/N) - cos(2*pi*k/N + pi/N)));
elseif (M == 3)
x = log2(abs(sin(pi*k/N).*(2*cos(2*pi/N) + 1) - sin(3*pi*k/N)));
elseif (rem(M, 4) == 0)
half = get_log2_zeros_multiplicative(N, M/2);
x = circshift(half, M/4) + circshift(half, -M/4);
elseif (rem(M, 4) == 1)
half = get_log2_zeros_multiplicative(N, (M-1)/2);
x = circshift(half, (M-1)/4+1) + circshift(half, -(M-1)/4) + get_log2_zeros_multiplicative(N, 1);
elseif (rem(M, 4) == 2)
half = get_log2_zeros_multiplicative(N, M/2);
x = circshift(half, (M-2)/4) + circshift(half, -((M-2)/4+1));
else
half = get_log2_zeros_multiplicative(N, (M-1)/2);
x = circshift(half, (M-3)/4+1) + circshift(half, -((M-3)/4+1)) + get_log2_zeros_multiplicative(N, 1);
endif
retval = x - round(max(x));
endfunction

function retval = geteig_multiplicative(N)
M = N-floor((N+2)/4)*2-1;
x = get_log2_zeros_multiplicative(N, M);
if (rem(M, 2) == 1)
x += get_log2_zeros_multiplicative(N, 1);
endif
x = ifftshift(2.^x);
firstzero = floor((N+10)/4);
x(firstzero:firstzero+M-1) = 0;
if (rem(N, 4) == 2 || rem(N, 4) == 3)
x0 = sum(x)/sqrt(N);
x1 = sum(x.*cos(2*pi*[0:N-1]/N))/sqrt(N);
s = x(1);
x .*= 2*cos(2*pi*[0:N-1]/N) + 2*(s - x1)/(x0 - s);
endif
retval = x/x(1);
endfunction

N = 512; x = geteig_multiplicative(N); plot([0:N-1], x - real(fft(x)/sqrt(N)), "x");

At $$N = 512$$, the error from the fast multiplicative algorithm is very low: Figure 6. The difference between the output of the fast multiplicative algorithm and its DFT, for $$N = 512$$.

For $$N = 65536$$, the absolute difference (not shown, analogous to Fig. 6) between the eigenvector calculated using the fast multiplicative algorithm and its DFT will be under 1e-12.

Hybrid multiplicative-iterative algorithm, $$O\big(\log(N)N\big)$$ time

The output of the multiplicative algorithm can still be refined using the iterative algorithm. In Octave:

function retval = geteig_hybrid(N)
if (N < 4)
retval = geteig_iterative(N, 0);
else
retval = geteig_iterative(N, 0, geteig_multiplicative(N));
endif
endfunction

N = 512; x = geteig_hybrid(N); plot([0:N-1], x - real(fft(x))/sqrt(N), "x"); Figure 7. The difference between the output of the fast hybrid multiplicative-iterative algorithm and its DFT, for $$N = 512$$, calculated in 9 refining iterations.

Using the iterative refinement, for $$N = 65536$$, the absolute difference (not shown, analogous to Fig. 7), between the eigenvector calculated using the fast multiplicative algorithm and its DFT, can be reduced from under 1e-12 to under 3e-16 in 23 iterations.

Conclusion

The multiplicative algorithm geteig_multiplicative(N) and the hybrid multiplicative-iterative algorithm geteig_hybrid(N) are very fast even for very large $$N$$, and are the recommended algorithms to use for generating these eigenvectors. The iterative refinement by the hybrid algorithm will be blind to any error due to Octave's FFT being done in double-precision floating point, which may in some applications give more error than the multiplicative algorithm which it uses to calculate the starting point of its iterations. The output of the multiplicative algorithm is guaranteed to be unimodal (Fig. 8), in the unwrapped sense, unlike the output of the hybrid algorithm. Figure 8. Base-10 logarithm of the unwrapped eigenvector calculated by the multiplicative algorithm (thick blue line) and its best-fit quadratic trendline (thin black line), for $$N = 65536$$. The logarithm of a Gaussian function is quadratic. The close match demonstrates that as $$N\to\infty$$, this type of an eigenvector approaches a Gaussian function. The first zeros would be in a plot like this at times $$\pm$$ floor((N+6)/4), for any $$N$$, in the current case at $$\pm 16385$$, meaning that both the Gaussian and the eigenvector take a dive to extremely small values already a long way away from the zeros.

• The pattern I have found is that my formula gives "eigenvectors in profile" for all N with odds alternating single sample peak/double sample peak, same with the evens. The single sample peaks can be rotated to DC and become true eigenvectors. The double sample peaks you can't. By adding them to a single sample shifted version of themselves, a single peak can be produced. Unfortunately (I was hoping), these are not true eigenvectors. RB-J's add the DFT trick for even functions works, but you don't get a minimal vector. I'm looking for a better theoretical foundation than Kuznetsov's. Aug 24 '20 at 11:53
• Yeah, as in my article, I need to do the shift and add for the even cases. This uses up a zero. At N/2 zero/nonzero split, the width of the DFT kernel is N/2+1, which means the problem is overdetermined and no soluton. After the add, they match and a solution is possible. Your iterations may converge faster starting with the double sample peaks added once over. The problem I am having is that DC centered is the way to go, but the other members of the $W_L$ family have peaks that are lower, so where to center? Most the efforts in literature concentrate on DC centered odd N. Aug 24 '20 at 13:01
• I'm still tracking, nice work. The N mod 4 = 2,3 cases are still the bugger for me as well. The extra roots are clearly the key, particularly your last formulation with 'a'. For your N=6 case, and N=15 it appears one is falling on $-\pi$ and $-\pi/2$ respectively, but that doesn't really make sense. Aug 25 '20 at 17:01
• @CedronDawg for $N=6$ the zero is at -3.130001029157937 and for $N=15$ it is at -1.555023702588202, give or take a few decimals, so not related to $\pi$ I guess. Aug 26 '20 at 5:58
• @CedronDawg symbolically for $N=6$ the convolution is by $[1, \sqrt{6}+1, 1]$ and for $N=10$ it is by $[1, \sqrt{10}/2 + \sqrt{5}/2 + \sqrt{2}/2 - 1/2, 1]$. (I did not look at $N=15$). The roots look even worse. Aug 26 '20 at 18:59