1
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Would it maybe be better to ask what is the closed-form expression of a windowed sinusoid? $$ x[n] = \cos(\omega_0 n + \phi) w[n] $$ where $w[n]$ is well defined. $\endgroup$ May 22, 2023 at 22:15
  • $\begingroup$ Well, depends what's "better". Sine is the ideal to show a full derivation, particularly for ones learning, also serves as a flexible reference. Windowing's fair as its own topic. $\endgroup$ May 22, 2023 at 23:46
  • 1
    $\begingroup$ I mean, you're applying a rectangular window anyway. So there will be sidelobes due to leakage. May as well make it a general window. $\endgroup$ May 23, 2023 at 0:12
  • $\begingroup$ Have you read the answer? The whole point is the rectangular window. A general window is a tremendous conceptual complication. $\endgroup$ May 23, 2023 at 0:14
  • $\begingroup$ I haven't gone through the answer. It looks like you're first dealing with DTFT and then, I presume, you will uniformly sample the DTFT to get the DFT. $\endgroup$ May 23, 2023 at 0:16

2 Answers 2

4
$\begingroup$

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 ) )
$\endgroup$
7
  • 1
    $\begingroup$ @OverLordGoldDragon: Reading the title of the question, there's no sampling or $t$ as a continuous time variable involved. I don't see the added value if we're just interested in an expression for the DFT of a discrete-time sinusoid. $\endgroup$
    – Matt L.
    May 24, 2023 at 16:53
  • $\begingroup$ @OverLordGoldDragon: Sampling rate is just relevant in relation to frequency, and this ratio is expressed by $\omega_0=2\pi f/f_s$ (if we even want to think of $x[n]$ as a sampled continuous-time sinusoid). $\endgroup$
    – Matt L.
    May 24, 2023 at 17:00
  • 2
    $\begingroup$ Concerning windows, I think we can't expect a closed-form solution for all types of windows. However, there could be one for windows defined using the cosine function, such as Hamming or Hanning. But anyway, the original question wasn't about windows, and you even said yourself "A general window is a tremendous conceptual complication". $\endgroup$
    – Matt L.
    May 28, 2023 at 11:03
  • $\begingroup$ Concerning vandalism, I don't take the time to check edits of other people's posts by whoever. We have the option of rolling back undesired edits, and I've done it myself a couple of times. As for the -2, I agree that downvoters should leave comments, but as you know it's not obligatory. If you're happy with your answer, why care about downvotes? $\endgroup$
    – Matt L.
    May 28, 2023 at 11:06
  • $\begingroup$ Don't give too much weight to votes. They're subjective and it's ok if some people don't like your answer (just 2 out of 137 who have seen this question). If it makes you feel better, take a look at this question. The accepted answer has 4 upvotes (and I think it's wrong), whereas my answer (which I think is correct) has zero votes. It's ok. $\endgroup$
    – Matt L.
    May 28, 2023 at 11:18
1
$\begingroup$

$$ \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,

  1. 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.
  2. Unwindowed STFT of sine, closed form solution (sliding FFT) - follows straightforwardly, using previous bullet's logic (yet even easier math), here.
  3. Energy behavior, $N/4$-symmetry, $1/x$ decay, asymmetry near DC & Nyquist - see "The modulus" (complex is also symmetric).
  4. 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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.