$x_a(t) = \cos(2\pi f_a t)$ was sampled with sampling period $T_s$. Plot the { spectrum | $N$-point DFT } of $x[n]$ ($f_a$, $T_s$ or $f_s$ given, $N$ given - whole number of periods or not).

Anyone can help?

  • 2
    $\begingroup$ Hi! plot with Matlab / Octave ? $\endgroup$ – Fat32 Nov 4 '18 at 21:41
  • $\begingroup$ @Fat32 It's enough to plot it just on a piece of paper. $\endgroup$ – emil Nov 4 '18 at 22:05
  • $\begingroup$ ok. So where is the difficult part ? $\endgroup$ – Fat32 Nov 4 '18 at 22:06
  • $\begingroup$ @Fat32 I don't know how to even start or how it should look like. $\endgroup$ – emil Nov 4 '18 at 22:07
  • 1
    $\begingroup$ ok ;-) the lab is tomorrow at least I believe ? you have plenty of time to learn the spectrum of $x[n]= \delta[n]$ ? I will prepare an answer though. $\endgroup$ – Fat32 Nov 4 '18 at 22:14

Your signal $x(t) = \cos(2 \pi f_a t)$ is an ideally defined infinite length signal, and the process of representing its DFT spectrum can be accomodated by the following windowing argument:

Let $x_d[n]=x(nT)=\cos(2 \pi f_a T_s n)$ be an ideally defined infinite length discrete-time sequence resulting of sampling ideally defined infinitely long continuous-time signal $x(t)=\cos(2\pi f_a t)$. And let a practically defined finite length truncated verison of $x_d[n]$ be

$$x[n] = x_d[n] w[n] $$

where $w[n]$ is a rectangular window of length $N$ samples and is used to truncate $x_d[n]$ according to:

$$ w[n] = \begin{cases} 1 &, 0 \leq n \leq N-1 \\ 0 &, \text{ otherwise } \\ \end{cases} $$

At this very point one can simply write down the $N$-point DFT sum equation as

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{ -j\frac{2\pi}{N} k n } = \sum_{n=0}^{N-1} \cos(2\pi f_a T_s n) e^{ -j \frac{2\pi}{N} k n } $$ and with the help of Euler identity for the cosine you get :

$$ \begin{align} X[k] &= \sum_{n=0}^{N-1} 0.5 \left( e^{j \frac{2\pi f_a} {F_s}n} +e^{-j \frac{2\pi f_a} {F_s}n} \right) e^{ -j\frac{2\pi}{N} k n } \\ & = 0.5 \sum_{n=0}^{N-1} e^{j \frac{2\pi f_a} {F_s}n} e^{ -j\frac{2\pi}{N} k n } + 0.5 \sum_{n=0}^{N-1} e^{-j \frac{2\pi f_a} {F_s}n} e^{ -j\frac{2\pi}{N} k n }\\ &= 0.5 \sum_{n=0}^{N-1} e^{j 2\pi ( \frac{f_a}{F_s}- \frac{k}{N}) n} + 0.5 \sum_{n=0}^{N-1} e^{-j 2\pi ( \frac{f_a}{F_s}+ \frac{k}{N}) n}\\ \end{align}$$

Now those exponential sums can be cast into a geometric sum formula and yield the following:

$$ \boxed{ X[k] = 0.5 \frac{ 1 - e^{j 2\pi ( \frac{f_a}{F_s}- \frac{k}{N}) N}}{1 - e^{j 2\pi ( \frac{f_a}{F_s}- \frac{k}{N}) } } + 0.5 \frac{ 1 - e^{-j 2\pi ( \frac{f_a}{F_s}+ \frac{k}{N}) N}} {1 - e^{-j 2\pi ( \frac{f_a}{F_s}+ \frac{k}{N}) } } ~~~, k = 0,1,2,...,N-1 }$$

Now depending on the relation between $f_a$, $F_s$ and $N$, you will typically have two distinct looking kind of spectrums.

First case assumes that the factor $ (\frac{f_a}{F_s} - \frac{k}{N})N = \frac{f_a N}{F_s} - k $ is an integer; equivalently $\frac{f_a N }{F_s}$ is an integer. When this is the case The spectrum will be full of zeros for all $k$ except at two points where $ k = \pm \frac{ f_s N }{ F_s} $. And at those points the magnitude will be $X[k] = N/2$.

For all other cases for which the ratio $\frac{f_a N}{F_s}$ is not an integer, then you will have a sinc type DFT spectrum and which is more typically observed in general. You should evaluate that for each $k$. Indees it's the DFT of a rectangular window (all ones), shifted to cosine frequency locations and sampled at those $k$ index frequencies.

A simple matlab /octave program confirms the results:

clc; clear all; close all;

fa = 2300;
Fs = 8000;
N = 81;

Ts = 1/Fs;
n = 0:N-1;

x = cos(2*pi*fa*n*Ts);      % Get truncated discrete-time sequence x[n].

X = fft(x,N);               % Get N-point DFT of x[n] via FFT function.

% Check our hand derivation:
k = 0:N-1 ;
Xk = 0.5* (  ( 1 - exp(-j*2*pi*( fa/Fs - k/N )*N))./(1 - exp(-j*2*pi*( fa/Fs - k/N ))))...
    + 0.5* (  ( 1 - exp(j*2*pi*( fa/Fs + k/N )*N))./(1 - exp(j*2*pi*( fa/Fs + k/N )))) ;

Xk(isnan(Xk)) = N/2;        % recover divide-by-zero when fa*N/Fs is an integer.

hold on
title([num2str(N),'-point DFT magnitude, |X[k]|, when the ratio (f_a N/F_s) is NOT an Integer']);
legend('via FFT','','by Hand');

with the following plots for $N=80$ and $N=81$ :

enter image description here

enter image description here

  • 1
    $\begingroup$ Check out my blog article at dsprelated.com/showarticle/771.php for how to greatly simplify your "by hand" formula. The result is equation (25). I have not seen this simplification anywhere else. The greatest advantage of the simplified formula is that it can be inverted. That is, the exact frequency of a single pure real tone in a DFT can be found from three bin values. I have been told repeatedly that an exact frequency formula is not possible, yet the first one I found is in my next article. Subsequent articles have other exact frequency formulas. $\endgroup$ – Cedron Dawg Nov 4 '18 at 23:34
  • $\begingroup$ @CedronDawg the three bin exact frequency computaion of a pure tone from its noise free DFT is also very nice example. Well appreciated your effort. Eventhough I already have it in my Matlab code base. You know that it requires very high SNR to be applied. And only single tone or you should modify the formula a bit ;-)... $\endgroup$ – Fat32 Nov 5 '18 at 8:51
  • $\begingroup$ I've added an answer to this post to show you what I mean by it being simplified. I am very interested in your claim that you have my formula in your MATLAB toolbox. Can you provide an online link to some documentation? The formula is actually quite robust in the presence of noise. Check out the reference to an independent evaluation made by Julien Arzi in the "Afterword". Jacobsen's Estimator is an approximation of my formula and has similar robustness. I have an even more robust formula in the pipeline. $\endgroup$ – Cedron Dawg Nov 5 '18 at 16:05
  • $\begingroup$ @CedronDawg Nooo I'm sorry, I didn't mean I had your formula in one of my matlab toolboxes, I wanted to say I have already coded a three DFT point exact sine frequency computation program by myself into my own codebase m-file (not in a toolbox) :-). Now reading your comment on "robust to noise" so probably they are not the same. I should read your formula better then... Becuse my simple derivation assumes that all DFT value are belonging to the pure sine signal , if there is noise, then DFT values will change hence the frequency estimation. How robust is your algorithm? SNR ? $\endgroup$ – Fat32 Nov 5 '18 at 16:36
  • $\begingroup$ The comparison paper can be found in the first link under the "Resources:" header in the gray box on the upper right on this page: tsdconseil.fr/log/scriptscilab/festim/index-en.html. You will see that mine was best in both low and high noise tests, especially the low noise cases. He provides a test case framework in SCILAB so you should be able to test yours in comparison. $\endgroup$ – Cedron Dawg Nov 5 '18 at 18:01

Let's start with the signal defined as a discrete sequence provided in Fat32's answer.

$$ x[n] = \cos( 2\pi \cdot f_a \cdot T_s n ) $$

With a simple variable substitution, this definition can be greatly simplified.

$$ \alpha = 2\pi \cdot f_a \cdot T_s $$

Here are the the units for the values in the equation:

$$ \frac{radians}{sample} = \frac{radians}{cycle} \cdot \frac{cycle}{second} \cdot \frac{seconds}{sample} $$

I find knowing the units tends to make equations a lot more understandable. The signal now looks like this:

$$ x[n] = \cos( \alpha n ) $$

Note that the units in the argument are really simple.

$$ \frac{radians}{sample} \cdot samples = radians $$

Next, let's precalculate two real values:

$$ U = x[N] - x[0] $$ $$ V = x[N-1] - x[-1] $$

When there are a whole number of cycles within the sample frame, both these values will be zero.

Lastly, define the relative location of each bin around the complex unit circle:

$$ \beta_k = k \cdot \frac{ 2\pi }{ N } $$

Fat32's "by hand" equation is now equivalent to this:

$$ X[k] = 0.5 \left[ \frac{ Ue^{ i\beta_k } - V }{ cos( \alpha ) - cos( \beta_k )} \right], ~~~ k = 0,1,2,...,N-1 $$

As you can plainly see, not only does this form require significantly fewer calculations to evaluate, but it also comprehensible qualitatively. The only complex valued term is in the exponent in the numerator.

The derivation can be found in my blog article DFT Bin Value Formulas for Pure Real Tones. In the article a variable for the amplitude of the signal and a phase value is also included. I also use a $1/N$ normalized DFT. I highly recommend you read the "Qualitative Analysis" section to see how the equation explains the well known behavior of a single pure tone in a DFT.

In the whole cycle case, the equation becomes $0/0$ at bins $k_f$ and $(N-k_f)$, where $k_f$ is the bin corresponding to the cycles per frame of the signal. Fat32's answer and my blog article explain how to deal with it.

I can pretty much guarantee that your professor has never seen this equation before.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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