Is the discrete Gaussian kernel an eigenfunction of the DFT?

Question

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?

Answer 1

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.

Answer 2

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

• The Zeroing Sine Family of Window Functions

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] &= 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[0] = 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

Answer 3

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. It is similar to Projection Onto Convex Sets (POCS) algorithms. 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[0] + x[1] = X[0] + a\,X[0] + X[0]$. 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.

Answer 4

So @AndreasH. said this:

... 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 ...

I think that's the answer and I'm gonna try to develop that into a constructive proof.

I'm starting with these definitions:

Fourier Transform: $$ \mathscr{F}\Big\{x(t)\Big\} \triangleq X(f) = \int\limits_{-\infty}^{\infty}x(t)\, e^{-j 2 \pi f t} \,\mathrm{d}t$$

and inverse: $$ \mathscr{F}^{-1}\Big\{X(f)\Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty}X(f)\, e^{j 2 \pi f t} \,\mathrm{d}f$$

Sampling function: $$ T \, \mathbf{III}_T(t) \triangleq T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) $$

We can uniformly sample $x(t)$ at the sample rate, $f_\text{s} \triangleq \frac{1}{T}$.

$$\begin{align} x_\text{s}(t) &= x(t) \cdot T \, \mathbf{III}_T(t) \\ &= x(t) \cdot T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} x(t) \, \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} x(nT) \, \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} x[n] \, \delta(t - nT) \\ \end{align}$$

Where $x[n] \triangleq x(nT)$.

It is also true that the sampling function is periodic and has a Fourier series.

$$\begin{align} T \, \mathbf{III}_T(t) &\triangleq T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\ &= \sum\limits_{k=-\infty}^{\infty} e^{j 2 \pi k f_\text{s} t} \\ \end{align}$$

Turns out that all of the Fourier series coefficients are 1. This means that the uniform sampled function is

$$\begin{align} x_\text{s}(t) &= x(t) \cdot T \, \mathbf{III}_T(t) \\ &= x(t) \cdot T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\ &= x(t) \sum\limits_{k=-\infty}^{\infty} e^{j 2 \pi k f_\text{s} t} \\ &= \sum\limits_{k=-\infty}^{\infty} x(t) \, e^{j 2 \pi k f_\text{s} t} \\ \end{align}$$

Accordingly, taking the continuous Fourier Transform, the spectrum of the sampled signal is

$$\begin{align} X_\text{s}(f) & \triangleq \mathscr{F} \Big\{ x_\text{s}(t) \Big\} \\ &= \mathscr{F} \left\{ \sum\limits_{k=-\infty}^{\infty} x(t) \, e^{j 2 \pi k f_\text{s} t} \right\} \\ &= \sum\limits_{k=-\infty}^{\infty} \mathscr{F} \Big\{ x(t) \, e^{j 2 \pi k f_\text{s} t} \Big\} \\ &= \sum\limits_{k=-\infty}^{\infty} X(f - k f_\text{s}) \\ \end{align}$$

And because of the Duality Theorem of the Fourier Transform, we know that the inverse applies if we uniformly sample in the frequency domain.

We also know that the Fourier Transform of the Gaussian is a Gaussian:

$$ \mathscr{F}\Big\{e^{-\pi t^2}\Big\} = e^{-\pi f^2} $$

$$ \mathscr{F}\Big\{e^{-\pi \alpha t^2}\Big\} = \frac{1}{\sqrt{\alpha}} e^{-\pi f^2/\alpha} $$

So now define $X_\text{s}(f)$ to be:

$$ X_\text{s}(f) \triangleq \sum\limits_{m=-\infty}^{\infty} \frac{1}{\sqrt{\alpha}} e^{-\pi (f-mf_\text{s})^2/\alpha} $$

That corresponds to this sampled Gaussian:

$$\begin{align} x_\text{s}(t) &= e^{-\pi \alpha t^2} \cdot T \, \mathbf{III}_T(t) \\ &= e^{-\pi \alpha t^2} \cdot T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} e^{-\pi \alpha (nT)^2} \, \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} x[n] \, \delta(t - nT) \\ \end{align}$$

Where $x[n] = e^{-\pi \alpha (nT)^2}$.

Now, so far $x_\text{s}(t)$ is continuous-time, it's sampled but still defined in continuous-time. And it's not (yet) periodic.

So now what happens if we sample the continuous-frequency and periodic function $X_\text{s}(f)$? Since it's periodic with period $f_\text{s}$, we have to sample it more densely with sampling interval of $\frac{1}{N}f_\text{s}$ in the $f$ domain.

$N \in \mathbb{Z} > 0$ is our standard DFT length and is a positive integer.

So using duality only changes every occurance of $-j$ with $+j$. Both have equal claim to squaring to be $-1$. Sampling $X_\text{s}(f)$ with sampling interval of $\frac{1}{N}f_\text{s}$ is:

$$\begin{align} \mathscr{F}^{-1} \Big\{ X_\text{s}(f) \tfrac{1}{N}f_\text{s}\sum\limits_{k=-\infty}^{\infty} \delta(f - \tfrac{k}{N}f_\text{s}) \Big\} &= \sum\limits_{m=-\infty}^{\infty} x_\text{s}(t-mNT) \\ &= \sum\limits_{m=-\infty}^{\infty} x(t-mNT) \cdot T \sum\limits_{n=-\infty}^{\infty} \delta(t - mNT - nT) \\ &= \sum\limits_{m=-\infty}^{\infty} x(t-mNT) \cdot T \sum\limits_{n=-\infty}^{\infty} \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} \sum\limits_{m=-\infty}^{\infty} x(t-mNT) \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} \sum\limits_{m=-\infty}^{\infty} x(nT-mNT) \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} \sum\limits_{m=-\infty}^{\infty} x[n-mN] \delta(t - nT) \\ &= T \sum\limits_{n=-\infty}^{\infty} \tilde{x}[n] \delta(t - nT) \\ \end{align}$$

$$\begin{align} \mathscr{F}^{-1} \Big\{ X_\text{s}(f) \tfrac{1}{N}f_\text{s}\sum\limits_{k=-\infty}^{\infty} \delta(f - \tfrac{k}{N}f_\text{s}) \Big\} &= \mathscr{F}^{-1} \Big\{ \tfrac{1}{N}f_\text{s}\sum\limits_{k=-\infty}^{\infty} X_\text{s}(f) \delta(f - \tfrac{k}{N}f_\text{s}) \Big\} \\ &= \mathscr{F}^{-1} \Big\{ \tfrac{1}{N}f_\text{s}\sum\limits_{k=-\infty}^{\infty} X_\text{s}(\tfrac{k}{N}f_\text{s}) \delta(f - \tfrac{k}{N}f_\text{s}) \Big\} \\ \end{align}$$