Unwindowed STFT of sine, closed form solution and insights (sliding FFT)

Question

I seek to calculate, mathematically, the unwindowed Short-Time Fourier Transform of

$$ \cos(2\pi f t + \phi) $$

i.e. any arbitrary real-valued sine: any frequency, duration, phase shift, and number of samples. The expression should reflect practical use, in that it's how we'd sample a sine in code. In effect, something like (MATLAB)

t = [0:N-1] * sampling_period
x = cos(2*pi * f * t + phi)
closed_form(N, f, phi, sampling_period, nfft, hop_size) == stft(x, nfft, hop_size) 

Answers can merely specify the steps required to perform this calculation, but it should be enough such that, the guided solution, expressed in code, validates against stft(x).

Once the solution is obtained, can it be used to reveal something new, or confirm known facts (e.g. symmetry, decay behavior)? Optionally, the result should be code-validated.

Answer 1

$$ \boxed{ \begin{align} & \texttt{STFT}_{M, H}\{\cos(2\pi f t + \phi)\}_{f(M/N)\notin\mathbb{Z}}[k, \tau] = \\ &\qquad\qquad\qquad \sin(\pi f_M) \left( \frac{U_\tau e^{j2\pi k/M} - V_\tau} {\cos(2\pi k/M) - \cos(2\pi f_M/M)} \right) \\ & U_\tau = \sin(\pi f_M + \phi_\tau),\\ & V_\tau = \sin(\pi f_M + \phi_\tau - 2\pi f_M/M),\\ & \phi_\tau = \phi + 2\pi (f_M/M)H\tau,\ f_M = f(M/N). \\ & \qquad H = \text{hop size},\ M = \text{segment size} = \texttt{nfft},\ N = \text{signal size} \\ & \qquad \texttt{STFT}_{M, H}\{x(t)\}[k, \tau] = \texttt{DFT}\{x(t + \tau H/N)_{:M}\}[k] \\ & \qquad x(s)_{:M} = x(s_{:M}),\ s_{:M} = [s[0], s[1], ..., s[M - 1]] \\ & \qquad \tau, H\in \mathbb{Z}; \ 0 \leq \tau \leq \lfloor{(N - M)/H\rfloor} + 1; \ 1 \leq H \leq N;\ M \leq N \\ & \qquad t = (1/N)[0, 1, ..., N - 1] \end{align} } $$

and (with same $t$ & others),

$$ \boxed{ \begin{align} & \texttt{STFT}_{M, H}\{\cos(2\pi f t + \phi)\}_{f(M/N)\in\mathbb{Z}}[k, \tau] = \\ & \qquad \frac{M}{2}\left(e^{j\phi_\tau} \delta [(k - f_M)_M] + e^{-j\phi_\tau}\delta [(k + f_M)_M]\right) \\ & (z)_M = z\ \text{mod}\ M \end{align} } $$

This implements 'valid'-mode (unpadded / don't hop outside of x) STFT, or also (with $N$ modified) STFT of a proper (e.g. 'reflect'-) padding of a sine that is within one sample of whole number of cycles.

Except for being precise on discrete vs continuous units, the result follows trivially from DFT of sine:

1. We just seek $\texttt{DFT}\{\cos(2\pi f(t + \tau) + \phi)\}$. Equating to $\texttt{DFT}\{\cos(2\pi f t + \phi_\tau)\}$ and solving for $\phi_\tau$ yields $\phi_\tau = \phi + 2\pi f\tau$. Instead of just $\tau$, we need $\tau/N$ so that $\tau$ is integer ($/N$ per how $t$'s defined) and hence refers to shift in samples and is a valid array index. Finally, accounting for hop brings us to $\tau H/N$. Replace $\phi$ with $\phi_\tau$ in the sine solution.
2. Recall that shorter DFT "sees" lower $f$ (see "Understanding SR vs Duration" here). So, replace $f$ with $f_M = f(M/N)$ in sine solution.
3. Replace $N$ with $M$ in sine solution (we're taking $M$-point DFTs)
4. Per 2, DFT's "integer $f$" is when $f(M/N)$ is integer, reproducing integer-$f$ sine DFT.

All other stuff is to exactly match STFT as we'd compute it:

• $H$, in math sense, is same as plugging in different $\tau$.
• Bounds on $\tau$ and $M$ make it 'valid', and on $H$ make it invertible.
• $s_{:M}$ stuff is a strictly correct way of saying "$M$-point DFT" without leaving $t$'s bounds.

For clarity: recall what $t$ is, and let $\tau = 1$, $H=2$, so now we're looking at the cosine at $[2, 3, ..., 2 + (M - 1)]$, so $\tau=0$ was fft(x[:M]) and $\tau=1$ is fft(x[2:M+2]).

In code (& $t$ offset note)

With sine_dft, this becomes the full code (np = numpy):

def sine_stft(N, M, H, f, phi):
    assert M <= N and 1 <= H <= N, (N, M, H)
    n_hops = (N - M)//H + 1
    out = np.zeros((M, n_hops), dtype='complex128')

    for tau in range(n_hops):
        phi_tau = phi + 2*np.pi*f*tau*H/N
        out[:, tau] = sine_dft(M, f*M/N, phi_tau)
    return out

For $t = [t_0, ...]$, that's sine_stft(N, M, H, f, phi + 2*pi*f*t[0]) (see "Effects on parameters" here).

Spectrogram

$$ \boxed{ \begin{align} & \left|\texttt{STFT}_{M, H}\{\cos(2\pi f t + \phi)\}_{f_M\in\mathbb{Z}, f_M\notin[0, M/2]}[k, \tau]\right| = \frac{M}{2}\big(\delta[(k - f_M)_M] + \delta[(k + f_M)_M]\big) \\ & \left|\texttt{STFT}_{M, H}\{\cos(2\pi f t + \phi)\}_{f_M\notin\mathbb{Z}}[k, \tau]\right| = \left|\sin(\pi f_M)\right| \frac{\sqrt{U_\tau^2 + V_\tau^2 - 2U_\tau V_\tau\cos(2\pi k/M)}} {\left|\cos(2\pi k/M) - \cos(2\pi f_M/M)\right|} \\ & U_\tau = \sin(\pi f_M + \phi_\tau),\\ & V_\tau = \sin(\pi f_M + \phi_\tau - 2\pi f_M/M),\\ & \phi_\tau = \phi + 2\pi (f_M/M)H\tau,\ f_M = f(M/N). \\ \end{align} } $$

and if we don't exclude the two edge cases,

$$ \boxed{ \begin{align} & \left|\texttt{STFT}_{M, H}\{\cos(2\pi f t + \phi)\}_{f_M\in\mathbb{Z}}[k, \tau]\right| = \\ &\qquad \frac{M}{2}\sqrt{ \delta[(k - f_M)_M] + \delta[(k + f_M)_M] + \delta[(k - f_M)_M] \delta[(k + f_M)_M](2\cos(2\phi_\tau)) } \\ \end{align} } $$

This is best interpreted via familiarity with sine DFT modulus, explored in-depth here.

(Note, the $U, V$ formulation is fully equivalent to taking modulus of the expression on top of this answer, explained in "Addendum: Original version vs Modified" here.)

Interpretation

STFT slides a window over input and takes its DFT. Here there's no window, so we're sliding DFT over a sine. This is exactly the same as sliding the sine instead and keeping DFT in same place. The operations are identical, but the interpretations aren't:

• "sliding DFT" is well-understood as, we're inspecting what the input is like over a certain interval, and using the spectrum to describe it - over various such intervals
• "sliding sine" is well-understood as, we're tracking how the spectrum of the sine changes - how it compares with the unshifted sine, what effect shifting is having. This isn't as relevant for the general input - a sine is stationary: fixed frequency and amplitude over time, so all that's changing is what portion of it we're looking at with a finite frame.

As a sine is being shifted in time by $\tau$,

• All imaginary bins are modulated by a sine of frequency $f$ as a function of $\tau$, with modulation amplitude and offset that depends on the sine.
• A given real bin modulated by a sine of frequency $f$ as a function of $\tau$, with modulation amplitude and offset that depends on the sine and the bin.

where "offset" = fixed phase, and "the sine" = original, unshifted $x(t)$. This is tailored to the second perspective, with an important application - see "Application: High SNR contamination" in Sine DFT article. An interpretation per first perspective can also be made, using time-frequency analysis - left to reader.

Insights

1. Unwindowed = Ideal for $f(M/N)\in\mathbb{Z}$

The formula doesn't lie, try N, M, f = 128, 32, 20 with any H & others.

2. Unwindowed = Ideal, for parameter retrieval (noiseless)

Since the inverse problem has been solved (see "Application: exact ..." in original article).

3. Beats near DC & Nyquist

$f$ is swept logarithmically over $[0.1, N/4]$, then backwards up to $N/2$, with $N, M = 512, 64$ and $\phi=0$.

The beats can be explained by sine DFT modulus (Ctrl + F "heavily phase-dependent"), where it was found that energy has the greatest $\phi$-dependence for $f$ near DC & Nyquist. Note, the GIF shows sweeping of $f$; sweeping of $\phi$ is the STFT itself (its time axis).

Directly reasoning from sine DFT alone doesn't generalize, however. The general case (any windowing, multi-component and AM-FM signals) is proven with time-frequency analysis: Why are there beats in spectrogram of sines?. Indeed, this beating persists for any windowing, and given STFT's time vs frequency symmetries - i.e. involving DFTs along frequency and time (explanation) - studying sine's STFT can reveal things about sine's DFT (which I've not done).

Full animation (epilepsy warning)

4. Sine Sliding FFT = ratio of frame edge-modulated & -offset roots of unity

An insight from the DFT sine solution ("Signal and Roots of Unity Formulation"):

$$ X[k] = \frac{(x[N] - x[0])R_k - (x[N - 1] - x[-1])}{(x[-1] + x[1])/x[0] - (R_k + R_{-k})} $$

where $R_k = e^{j2\pi k/N}$, and indexing outside of $[0, N-1]$ is of sine's continuation, not circular.

The DFT knows no such thing as "continuation" beyond what it's given, but STFT's frames do. This can interpret, and we can observe, as follows:

1. Guaranteed imaginary part: since $x[N] = x[0]$ can't persist for all $\tau$. Exceptions: "all $\tau$" is one $\tau$; $H=M$ and $f_M$ is integer (for certain $\phi$).
2. Rare $|\texttt{STFT}[:, \tau_0]| = \texttt{STFT}[:, \tau_0]$ (imag=0): since $x[N] = x[0]$ is possible. "$[:, \tau_0]$" means $\forall k$.
3. Difference of edges controls decay rate, since it controls the numerator (point 4 here).
4. Difference of edges controls decay behavior of peaks near DC & Nyquist (point 1 here)
5. Ratio of starting points controls peak location & intensity, since it controls the denominator. It's sum of $x[\pm 1]$ divided by their midpoint.
6. Difference of edges controls imaginary part (its strength relative to real part), since $x[N] - x[0]$ controls the only imag term, and the denominator affects real & imag same.
7. Fully determined by 5 points per frame, 3 at start and 2 at end, one out-of-frame per edge.

Citation

This work can be cited in one or two parts:

John Muradeli, 2023. Unwindowed STFT of sine, closed form solution and insights (sliding FFT). URL: https://dsp.stackexchange.com/a/88806/50076

Cedron Dawg, 2015. DFT Bin Value Formulas for Pure Real Tones. URL: https://www.dsprelated.com/showarticle/771.php

If must cite only one, cite first.

Code validation

Thoroughly validated against scipy.signal.stft, available on Github.

Note, I validated via sine_dft, so the exact LaTeX above wasn't coded and assumes I got all the substitutions right. Also, sine_dft follows the original sine solution, while what's presented in this article is modified; they're equivalent, see "Addendum" in original article. I was careful in handling these correctly, but if it matters, double-check.

Animation code exists, can't share at the moment.