DFT of a sine, closed form solution and insights

Question

I seek to calculate, mathematically, the Discrete Fourier Transform,

$$ \texttt{DFT}\{x\}[k] = \sum_{n=0}^{N - 1} x[n] e^{-j2\pi k n / N} $$

of 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) == fft(x) 

Answers are free to formulate their own expression for $x[n]$ and $x(t)$. An incomplete example is $x(t) = \cos(\omega t + \phi)$ (missing sampling rate, etc). Answers can also merely specify the steps required to perform this calculation, but it should be enough such that, the guided solution, expressed in code, validates against fft(x). Infinite sums don't qualify.

Once the solution is obtained, can it be used to reveal something new, or confirm known facts (e.g. symmetry, decay behavior)? This part is optional (including a yes/no). Also optionally, the result should be code-validated.

Answer 1

The length $N$ DFT of the discrete-time sinusoid

$$x[n]=\cos(\omega_0n+\phi),\qquad n=0,1,\ldots,N-1\tag{1}$$

is given by

$$X[k]=\sum_{n=0}^{N-1}\cos(\omega_0n+\phi)e^{-j2\pi nk/N},\qquad k=0,1,\ldots,N-1\tag{2}$$

Eq. $(2)$ can be rewritten as

$$\begin{align}X[k]&=\frac12\sum_{n=0}^{N-1}\left[e^{j(\omega_0n+\phi)}+e^{-j(\omega_0n+\phi)}\right]e^{-j2\pi nk/N}\end{align}\tag{3}$$

The sum over the first term inside brackets in $(3)$ can be computed as follows:

$$\begin{align}\sum_{n=0}^{N-1}e^{j(\omega_0n+\phi)}e^{-j2\pi nk/N}&=e^{j\phi}\sum_{n=0}^{N-1}e^{jn\left(\omega_0-\frac{2\pi k}{N}\right)}\\&=e^{j\phi}\frac{1-e^{jN\omega_0}}{1-e^{j\left(\omega_0-\frac{2\pi k}{N}\right)}}\\&=e^{j\phi}\frac{e^{j\frac{N\omega_0}{2}}}{e^{j\left(\frac{\omega_0}{2}-\frac{\pi k}{N}\right)}}\cdot\frac{\sin\left(\frac{N\omega_0}{2}\right)}{\sin\left(\frac{\omega_0}{2}-\frac{\pi k}{N}\right)}\\&=e^{j\left(\frac{N-1}{2}\omega_0+\frac{\pi k}{N}+\phi\right)}\frac{\sin\left(\frac{N\omega_0}{2}\right)}{\sin\left(\frac{\omega_0}{2}-\frac{\pi k}{N}\right)}\tag{4}\end{align}$$

Note that in $(4)$ we've assumed that $\omega_0\neq 2\pi m/N$, $m\in\mathbb{Z}$, because otherwise there would be a division by zero for $(k-m)_N=0$, where the $(\cdot)_N$ denotes a modulo $N$ operation. For all other values of $k$ the sum equals zero in that case. Hence, for $\omega_0 = 2\pi m/N$, $m\in\mathbb{Z}$, the top left sum in Eq. $(4)$ is given by

$$\sum_{n=0}^{N-1}e^{j(\omega_0n+\phi)}e^{-j2\pi nk/N}=Ne^{j\phi}\delta[(k-m)_N]\tag{5}$$

The sum over the second term in brackets in Eq. $(3)$ is obtained in a completely analogous manner, leading to the final result

$$X[k]=\frac12e^{j\left(\frac{N-1}{2}\omega_0+\frac{\pi k}{N}+\phi\right)}\frac{\sin\left(\frac{N\omega_0}{2}\right)}{\sin\left(\frac{\omega_0}{2}-\frac{\pi k}{N}\right)}+\frac12e^{-j\left(\frac{N-1}{2}\omega_0-\frac{\pi k}{N}+\phi\right)}\frac{\sin\left(\frac{N\omega_0}{2}\right)}{\sin\left(\frac{\omega_0}{2}+\frac{\pi k}{N}\right)}\tag{6}$$

which is valid for $\omega_0\neq2\pi m/N$, $m\in\mathbb{Z}$. For $\omega_0=2\pi m/N$, $m\in\mathbb{Z}$, the result is

$$X[k]=\frac{N}{2}e^{j\phi}\delta[(k-m)_N]+\frac{N}{2}e^{-j\phi}\delta[(k+m)_N]\tag{7}$$

This Matlab/Octave code can be used to verify Eqs $(6)$ and $(7)$:

% x[n] = cos( 2*pi*f0*n + phi ), n=0,1,...,N-1
% choose these 3 parameters
N = 256;
phi = .2;
f0 = 8.3/N; % (integer number)/N gives integer number of cycles inside window

w0 = 2*pi*f0;

% check for integer number of cycles inside window
intcyc = 0;
m = f0*N;
if( floor( m ) == m ), intcyc = 1; endif

% evaluate analytical formula
if ( intcyc ),
    X = zeros(N,1);
    if( mod( m, N ) == 0 ),
        X( 1 ) = N * cos(phi);
    else
        idx = mod( m, N ) + 1;
        X( idx ) = N/2 * exp( 1i*phi );
        idx = mod( -m, N ) + 1;
        X( idx ) = N/2 * exp( - 1i*phi );
    endif
else
    k = ( 0 : N-1 )';
    X = .5 * exp( 1i * ( (N-1)/2*w0 + pi*k / N + phi ) ) * sin( N*w0/2 ) ./ sin( w0/2 - pi*k/N );
    X = X + .5 * exp( - 1i * ( (N-1)/2*w0 - pi*k / N + phi ) ) * sin( N*w0/2 ) ./ sin( w0/2 + pi*k/N );
endif

% compare with FFT
n = ( 0 : N-1 )';
x = cos( w0*n + phi ); 
X2 = fft( x, N );

max( abs( X2 - X ) )

Answer 2

$$ \boxed{ \texttt{DFT}\{\cos(2\pi f t + \phi)\}_{f\notin\mathbb{Z}}[k] =\\ \qquad \sin(\pi f) \left( \frac{\sin(\pi f + \phi)e^{j2\pi k/N} - \sin(\pi f + \phi - 2\pi f/N)} {\cos(2\pi k/N) - \cos(2\pi f/N)} \right) } $$

where $t=\frac{1}{N}[0, 1, ..., N - 1]$ - and,

$$ \boxed{ \texttt{DFT}\{\cos(2\pi f t + \phi)\}_{f\in\mathbb{Z}}[k] = \frac{N}{2}\left(e^{j \phi} \delta [(k - f)_N] + e^{-j\phi}\delta [(k + f)_N]\right) } $$

where $(z)_N = z\ \text{mod}\ N$. Solution by Cedron Dawg, derivation included. Code validation at bottom.

Analysis utility

The expression is excellent for understanding the effects of input parameters on spectrum's behavior. It can be reduced to, and understood in terms of, four distinct terms:

$$ X[k] = \frac{1}{2} K \big[U e^{j2\pi k/N} - V\big] $$

where

$$ \begin{aligned} & K = \frac{1}{{\cos(2\pi f/N) - \cos(2\pi k/N)}} \\ & U = \cos(2\pi f + \phi) - \cos(\phi) \\ & V = \cos(2\pi f + \phi - 2\pi f/N) - \cos(\phi - 2\pi f/N) \end{aligned} $$

and $e^{j2\pi k/N}$ is the sole complex-valued term. One may notice, this doesn't match the top equation - in short, that's intended and they're equivalent (see Addendum). Parametrically, we observe:

• $K$ and $e^{j2\pi k/N}$ are the only $k$-dependent terms
• $f$ affects all terms
• $U, V$ have the mathematical form $\cos(A + B) - \cos(B)$.

Let "the denominator" mean $K$'s denominator. From these, we can infer (or confirm):

1. For integer $f$, $U=V=0$ and all bins are zero except $k=f$, where the denominator is also zero and yields indeterminate form.
2. As $k \rightarrow f$, $K$ is by far the most dominant term, hence why the DFT peaks for $k$ closest to $f$.
3. For $k$ deviating from $f$, the denominator grows and coefficients taper away from the peak.
4. When $f$ is close to an integer, $U$ and $V$ are small, and all bins are close to zero except for those near $f$, since there $K$ kicks in. $U, V$ are small since $\cos(A + B) = \cos(B)$ for $A$ that's a multiple of $2\pi$.
5. Near $f$, as $k$ goes from below $f$ to above it, the sign in the denominator changes, while the numerator stays roughly same, hence the nearly symmetric behavior around peaks.
6. For "more fractional" $f$ (far from being integer / closer to .5, 1.5, ...), all bins are larger (except peaks), since $U, V$ are larger per $\cos(A + B) \approx \cos(B)$ no longer holding ($A$ is far from integer multiple of $2\pi$).
7. For "more fractional" $f$, the peaks are smaller, since $K$ is now much smaller for $k \rightarrow f$ (note $1/x$-like behavior), overwhelming $U, V$.
8. $e^{j2\pi k/N}$ is the only complex-valued term, and it's also identical to the iDFT's kernel, which is well-known to conjugate above $N/2$. It clearly follows, $X[k] = X^{*}[N - k]$, and that (nearly) all observations for $f$ also hold for $N - f$. Besides "well-known", it's easily confirmed: $e^{j2\pi (N - k)/N} = e^{j2\pi (1 - k/N)} = e^{j2\pi}e^{-j2\pi k/N} = e^{-j2\pi k/N}$, thus $e^{j2\pi(N - k)/N} = (e^{j2\pi k/N})^*$.
9. Since $e^0 = 1$ and $e^{j2\pi (N/2)/N} = -1$, eliminating the only imaginary part, DC and Nyquist are real-valued.

and, with some work,

10. Spectrum varies sinusoidally with time-shifts, with frequency same as input's. Rewriting, $\cos(2\pi ft + \phi) = \cos(2\pi f(t + \tau))$, solving yields $\phi = 2\pi f \tau$. So, replace $\phi$ with $2\pi f \tau$ everywhere, and treat everything as constant except $\tau$; the only variables now are $U, V$. We have $U = \cos(2\pi f\tau + l_0) - \cos(2\pi f\tau)$ and $V = \cos(2\pi f\tau + l_1) - \cos(2\pi f\tau + l_2)$. Sum of sines of same frequency is another sine of same frequency, so $U = a\cos(2\pi f\tau + p)$ and $V = b\cos(2\pi f\tau + q)$. $V$ is real-valued, so it doesn't interact with the cisoid, so imaginary is $V(K/2)\sin(2\pi k/N)$. What we're really after is $\cos(2\pi f(t + \tau) + \phi)$, but this just changes $a, b, p, q$ - hence, $\boxed{\Im m\{X_\tau[k]\} = b\cos(2\pi f\tau + p)\Im m\{X[k]\}}$, with $x_\tau(t) = x(t + \tau)$. For the real part, it's same, except its version of $b, p$ are $k$-dependent. See "Proof: shifting" here.
11. Unwindowed STFT of sine, closed form solution (sliding FFT) - follows straightforwardly, using previous bullet's logic (yet even easier math), here.
12. Energy behavior, $N/4$-symmetry, $1/x$ decay, asymmetry near DC & Nyquist - see "The modulus" (complex is also symmetric).
13. Exact $f, \phi, A$ recovery from $X[k]$ - see below.

1-9 are also found in the solution article.

The modulus

Energy behavior, $N/4$-symmetry, $1/x$ decay, asymmetry near DC & Nyquist - and other insights, are revealed by studying $|X|$: DFT modulus of sine, closed form solution and insights.

Application: High SNR contamination

When testing an electronic measurement device, we can feed a near-ideal sine. Such sines can also be observed in other highly sensitive settings, like detecting gravitational waves (Einstein@Home), or working with other fixed-frequency radiation.

If we try to measure such SNR in a streamed manner, via scipy.signal.periodogram, we'll find one of the two:

The application is two-fold:

1. Understanding the problem: an engineer's time is valuable. This problem was asked right here on DSP.SE, and it took many collective hours, my own included, to figure out firstly why there's any shift-dependence at all, let alone why it's sine or sine-like. I eventually figured it out, but couldn't prove it - from the DFT solution, the answer, and the fix, follow easily. The problem was a flawed application of, or a maybe-flaw in, the scipy method. (Cause of problem follows easily, but proof of SNR behavior takes some work, done in "The modulus".)
2. Solving the problem, competitively: windowing performs excellent at reducing leakage. If it's not sufficient, we can enhance - the exact $f, \phi, A$ extraction (see below) performs phenomenally with such low noise, by which we can reproduce the noiseless spectrum, which not only is good by itself, but works together with windowing and other methods. ($\phi$ may be suboptimal, but a sweep yields energy upper bound.)

Application: Exact $f, \phi, A$ recovery

Traditional wisdom says, a finite sine is aliased, with DFT being the result of infinite overlap-addings of sinc, which are impossible to untangle. Traditional wisdom is wrong.

Not only is frequency recoverable exactly using the DFT, but so are amplitude ($A$) and phase. Oh, and it's doable with two bins: A Two Bin Solution, by Cedron Dawg. Using more bins, for retrieving frequency, it's also robust to noise - Exact Frequency Formula ..., by same author.

What's the buzz over, if unfamiliar? 1 and 2 have been open on Stack Exchange for 12 years, with multiple responders, and the best result is an approximation using STFT. Just two examples.

Note, it's not about being best at estimating, for which there's modern methods. It's about being exact with general $f, \phi, A$ when possible (noiseless), where it's thought to be "never" with DFT - forcing reevaluation of related theory and its implications.

Greater significance

Enhancing existing methods: the sine is a building block of all signals. Two examples were shown where a closed form solution is directly practically useful - it's just two examples. While I can think of others (some below), up to this point, the thinking and applying has been done almost exclusively by one person - the original author; hopefully this post invites further interest.

Ameliorating "spectral leakage": leakage is treated as a permanent corruption introduced by aliasing, for which there are only approximate remedies. Exact recovery of source sine parameters straight from the spectrum flips that on its head. While general signal leakage remains problematic, revisiting the assumptions in its treatment may produce superior treatments.

Machine Learning / Super Resolution: one source of imaginary part, one $k$-dependent in numerator and denominator. This computationally minimal form is highly extendable. It solves the inverse problem for $f$ with matrix manipulations of bin values - no nonlinearities. "Prior injection" can improve resolving multi-tone signals - retrieving closely-spaced $f_0, f_1, ...$. An example with simply DFT, using Compressive Sensing, is provided here (Royi Avital).

The aforementioned solution to $f$ can already estimate multi-tones with an iterative closest-fit algorithm, and is also enhanceable. Some methods can be enhanced to be noise robust - Cedron's one-$f$ solution is already robust, by exploiting the sine solution's mathematical form to reduce influence of irrelevant variability. This is key: a key factor in effectiveness of "enhancements" is how well the original method untangles variability. This is what much of "feature extraction" is about, and how it improves ML performance. A proper extension may claim State of the Art.

---

Interactive: Sine Solution (try live!)

(1/N-normalized DFT; x-axis is k/N)

Interactive: Any Signal (real-valued) (try live!)

Interactive: Any Signal (try live!)

---

Main reading is done. Code validation is at bottom.

---

Appendix A: Numerically precise version

Alt solution (same author):

$$ \boxed{ \texttt{DFT}\{\cos(2\pi f t + \phi)\}_{f\notin\mathbb{Z}}[k] = \\ \frac{1}{2} \frac{\sin(\pi f)}{\sin(\pi f/N - \pi k/N)} \left( \frac{\sin(\pi f + \phi)e^{j2\pi k/N} - \sin(\pi f + \phi - 2\pi f/N)} {\sin(\pi f/N + \pi k/N)} \right) } $$

This concerns $f \approx \text{integer}$. $a / (bc)$ is more precise than $a/(b - c)$, where $a, b, c \approx 0$. The original solution has the latter form: the difference of small numbers suffers more error than product of small numbers. More can be said, this is the gist.

Appendix B: Minimal code validation

Python:

from numpy import fft, allclose, pi, sin, cos, exp, arange

def sine_dft(N, f, phi):
    k = arange(N)
    num = sin(pi*f + phi)*exp(1j*2*pi * k/N) - sin(pi*f + phi - 2*pi*f/N)
    den = cos(2*pi*k/N) - cos(2*pi*f/N)
    return sin(pi*f) * num / den

N, f, phi = 131, 3.213, 1.2
x = cos(2*pi*f*arange(N)/N + phi)

assert allclose(sine_dft(N, f, phi), fft.fft(x))

MATLAB/Octave:

N = 131;
f = 3.213;
phi = 1.2;

k = (0:N-1);
num = sin(pi*f + phi)*exp(1j*2*pi * k/N) - sin(pi*f + phi - 2*pi*f/N);
den = cos(2*pi*k/N) - cos(2*pi*f/N);
sine_dft = sin(pi*f).* num ./ den;

x = cos(2*pi*f*(0:N-1)/N + phi);

assert(norm(sine_dft - fft(x)) < 1e-11)

Full code ("Code validation") also handles integer $f$.

Appendix C: Sampling Rate, Duration, $t$-offset

For physical/continuous interpretation and handling any $t$. Explained here. In short:

• Doubled duration: $f$ doubled
• Doubled sampling rate: $f$ unchanged
• Non-zero $t[0]$: $\phi \rightarrow \phi + 2\pi f t[0]$ (using old $f$)
• Duration, sampling rate: $T = (t[1] - t[0])\cdot N$, $S = 1 / (t[1] - t[0])$, $S = N/T$

In one equation: $\texttt{sine_dft}\{N, fT, \phi + 2\pi f t[0]\}$. In SR terms, $fT \rightarrow f(N/S)$.

Addendum: Original version vs Modified

$$ \boxed{ \texttt{DFT}\{\cos(2\pi f t + \phi)\}_{f \notin \mathbb{Z}}[k] = \frac{1}{2}\frac{1}{{\cos(2\pi f/N) - \cos(2\pi k/N)}} \left[ \\ \ \ (\cos(2\pi f + \phi) - \cos(\phi)) e^{j 2\pi k/N} - \left( \cos(2\pi f + \phi - 2\pi f/N) - \cos(\phi - 2\pi f/N) \right) \right] \\ } $$

is what $U, V$ expand to, and is the original solution.

In short, I found the other version after I already finished all the work in this post and ones I reference. Original's superior for some purposes, and modified for others. I elaborate here, under "Addendum: Original version vs Modified".

---

Citation

This work can be cited in one or two parts:

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

John Muradeli, 2023. DFT of a sine, closed form solution and insights. URL: https://dsp.stackexchange.com/a/88365/50076

Under "Analysis Utility", most points (i.e. not 10, 11, 12) repeat Cedron, so if citing for those or for "Application: Exact ...", only cite Cedron. The solution itself is also credited to Cedron. If citing for my insights, or something significant from other sections, please also cite this article.

Code validation

Available at Github.