Discrete Signal Processing (DSP):

Discrete Fourier Transformation (DFT)

Table of Contents





1. Eigenanalysis

1.1. Eigenvectors and Eigenvalues

Given a square matrix $A$, an eigenvector $v$ and corresponding eigenvalue $\lambda$ satisfy:


$$ Av = \lambda v $$

Extending to a full basis of eigenvectors, we can write:


$$ AV = V\Lambda $$

where $V = [v_0 \mid v_1 \mid \cdots \mid v_{N-1}]$ is the matrix of eigenvectors, and $\Lambda$ is the diagonal matrix of eigenvalues:


$$ \Lambda = \begin{bmatrix} \lambda_0 & & & \\ & \lambda_1 & & \\ & & \ddots & \\ & & & \lambda_{N-1} \end{bmatrix} $$

This leads to the eigen-decomposition of the matrix:


$$ A = V \Lambda V^{-1} \qquad \text{or equivalently} \qquad V^{-1} A V = \Lambda $$

This decomposition reveals how the matrix $A$ acts on a basis of its eigenvectors: it scales each one by a corresponding eigenvalue.



1.2. Eigenanalysis of LTI Systems (Finite-Length Signals)

For signals of length $N$, a linear time-invariant (LTI) system is represented by a circulant matrix $H \in \mathbb{C}^{N \times N}$ whose entries are determined by the impulse response $h[n]$:


$$ [H]_{n,m} = h[(n - m)_N] $$

Our goal is to find the eigenvectors and eigenvalues of such matrices.


Harmonic Sinusoids as Eigenvectors

The eigenvectors of any circulant matrix (and hence any finite-length LTI system) are known to be the harmonic sinusoids:


$$ s_k[n] = \frac{1}{\sqrt{N}} e^{j\frac{2\pi}{N}kn}, \quad \text{for } k = 0, \dots, N - 1 $$

To prove this, consider the circular convolution between $s_k[n]$ and $h[n]$:



$$\begin{align*} y[n] = s_k[n]\, \ast h[n] \, &= \, \sum \limits_{m=0}^{N-1} s_k[(n-m)_N]\, h[m]\, = \, \sum \limits_{m=0}^{N-1} \frac{e^{j\frac{2 \pi}{N}k(n-m)_N}}{\sqrt{N}}\, h[m]\\ \\ & = \, \sum \limits_{m=0}^{N-1}\frac{e^{j\frac{2 \pi}{N}k(n-m)}}{\sqrt{N}}\, h[m]\, =\, \sum \limits_{m=0}^{N-1}\frac{e^{j\frac{2 \pi}{N}kn}}{\sqrt{N}}e^{-j\frac{2 \pi}{N}km}\, \, h[m]\\ \\ & = \, \left(\sum \limits_{m=0}^{N-1}e^{-j\frac{2 \pi}{N}km}\, h[m]\,\right)\frac{e^{j\frac{2 \pi}{N}kn}}{\sqrt{N}} \, = \, \lambda _k \, s_k[n] \end{align*}$$

Thus, $s_k[n]$ is an eigenvector of $H$, and its corresponding eigenvalue is:


$$ \lambda_k = \sum_{n=0}^{N-1} h[n] e^{-j\frac{2\pi}{N}kn} = \langle h, s_k \rangle = H_u[k] $$

This eigenvalue $\lambda \in \mathbb{C}$ is interpreted as the frequency response of the system at frequency index $k$, measuring how strongly the system responds to the complex exponential $s_k$.



1.3. Eigenvector Matrix of Harmonic Sinusoids

We can stack the $N$ orthonormal eigenvectors $\{s_k[n]\}_{k=0}^{N-1}$ into a unitary matrix $S$:


$$ \begin{align*} S &= \left[s_0 \quad s_1 \quad \cdots \quad s_{N-1} \right] \\ \\ &= \begin{bmatrix} s_0[0] & s_1[0] & \cdots & s_{N-1}[0] \\ s_0[1] & s_1[1] & \cdots & s_{N-1}[1] \\ \vdots & \vdots & & \vdots \\ s_0[N-1] & s_1[N-1] & \cdots & s_{N-1}[N-1] \end{bmatrix} \end{align*} $$

Each column $s_k$ is defined by:


$$ s_k[n] = \frac{1}{\sqrt{N}} e^{j\frac{2\pi}{N}kn} $$

This matrix $S$ forms the discrete Fourier transform (DFT) basis, and diagonalizes any circulant matrix $H$:


$$ H = S \Lambda S^{-1} $$

This result is fundamental in signal processing, as it connects time-domain convolution with frequency-domain multiplication via the DFT.


The following example visualizes the complex matrix $S$ by separately illustrating its real and imaginary components.

In [ ]:
%plot -s 800,500

% visualize S matrix

N = 16;
n = 0:N-1;
k = 0:N-1;

S = 1/sqrt(N)*exp(1j*2*pi/N).^(n'*k);

subplot(1,2,1), imagesc(n,k,real(S)); colormap('jet'); axis('square');
title('real part'), xlabel('k'), ylabel('n')
subplot(1,2,2), imagesc(n,k,imag(S)); colormap('jet'); axis('square');
title('imaginary part'), xlabel('k'), ylabel('n')










    

































2. Signal Decomposition Using Harmonic Sinusoids
Signal decomposition refers to expressing a signal as a sum of fundamental components. To understand this idea, we return to linear algebra and explore vector spaces, basis vectors, and orthonormal transformations. Signal transforms, which map signals between different representations or "domains," are naturally described in terms of basis expansions.




2.1. Basis
A basis $\{b_k\}$ for a vector space $V$ is a set of vectors that satisfies:


    

  • Span: Every vector $x \in V$ can be written as a linear combination of the basis vectors:
    • 




$$
x = \sum_{k=0}^{N-1} \alpha_k b_k = \alpha_0 b_0 + \alpha_1 b_1 + \cdots + \alpha_{N-1} b_{N-1}
$$


    

  • Linear independence: No vector in the set can be written as a linear combination of the others.
    • 





We can organize the basis vectors as columns of a matrix:



$$
B = \left[b_0 \mid b_1 \mid \cdots \mid b_{N-1} \right]
$$


Using this matrix, a signal can be represented compactly as:



$$
\begin{align*}
x &  =\, \alpha_0 b_0 + \alpha_1 b_1 +\, \cdots\, + \alpha_{N-1} b_{N-1} \, =\,  \sum_{k=0}^{N-1} \alpha_k b_k \,\\\\
&=\, \left[b_0 \mid b_1 \mid b_2 \mid \cdots \mid b_{N-1} \right]\begin{bmatrix}
\alpha_0\\
\alpha_1\\
\vdots\\
\alpha_{N-1}
\end{bmatrix} = \, B\, \alpha
\end{align*}
$$


where $\alpha = [\alpha_0 \; \alpha_1 \; \cdots \; \alpha_{N-1}]^T$ is the coefficient vector.

















2.2. Orthonormal Basis
A set of vectors $\{b_k\}_{k=0}^{N-1}$ is an orthogonal basis for a vector space $V$ if the vectors are mutually orthogonal:



$$
\langle b_k, b_l \rangle = 0 \quad \text{for all } k \neq l
$$


If, in addition, each vector has unit norm:



$$
\| b_k \|_2 = 1 \quad \text{for all } k
$$


then $\{b_k\}$ is called an orthonormal basis.




In this case, the basis matrix $B = [b_0 \mid b_1 \mid \cdots \mid b_{N-1}]$ is unitary, meaning:



$$
B^H B = I \quad \Longrightarrow \quad B^{-1} = B^H
$$


This orthonormality condition expands to:



$$
\begin{align*}
B^H B &=
\begin{bmatrix}
b_0^H\\
b_1^H\\
\vdots\\
b_{N-1}^H
\end{bmatrix}
\begin{bmatrix}
b_0 & b_1 & \cdots & b_{N-1}
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & \cdots & 0 \\
0 & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 1
\end{bmatrix}
= I
\end{align*}
$$


This identity plays a central role in signal representation, enabling simple analysis and synthesis of signals using inner products and linear combinations.

















2.3. Signal Representation with Orthonormal Bases
Given an orthonormal basis matrix $B$, any signal $x \in \mathbb{C}^N$ can be represented as:



$$
\begin{align*}
x &= B \alpha = \sum_{k=0}^{N-1} \alpha_k b_k \quad &\text{(synthesis)} \\\\
\alpha &= B^H x \quad \text{or} \quad \alpha_k = \langle x, b_k \rangle \quad &\text{(analysis)}
\end{align*}
$$


    

  • Synthesis: The signal $x$ is constructed as a linear combination of the basis vectors $b_k$ weighted by coefficients $\alpha_k$.
    

    • 

  • Analysis: The coefficients $\alpha_k$ are obtained by projecting the signal $x$ onto the basis vectors, capturing how much each basis vector contributes to the signal.
    

    • 





This decomposition is powerful:


    

  • The synthesis equation reconstructs the signal using the basis.
    • 

  • The analysis equation extracts the signal’s coordinates in the basis.
    • 





Since $B$ is unitary ($B^H = B^{-1}$), the transformation between $x$ and $\alpha$ is perfectly invertible.

















2.4. Harmonic Sinusoids as an Orthonormal Basis
The complex exponential signals, or harmonic sinusoids, form an orthonormal basis for $\mathbb{C}^N$. Define the $k$-th basis vector as:



$$
s_k[n] = \frac{1}{\sqrt{N}} e^{j \frac{2\pi}{N} kn}, \qquad 0 \leq k < N, \quad 0 \leq n < N
$$


Stacking these vectors as columns gives the $N \times N$ matrix:



$$
S = \left[ s_0 \mid s_1 \mid \cdots \mid s_{N-1} \right]
$$


Then the matrix $S$ satisfies:



$$
S^H S = I \quad \Rightarrow \quad \langle s_k, s_l \rangle = \delta_{k,l}, \quad \|s_k\|_2 = 1
$$


This shows that the set $\{s_k\}$ is an orthonormal basis.




These sinusoids are fundamental in signal processing and form the basis of the Discrete Fourier Transform (DFT), enabling analysis of signals in the frequency domain.




















3. Discrete Fourier Transform (DFT)
Jean Baptiste Joseph Fourier proposed a powerful idea: any signal can be represented as a linear combination of sinusoids.




3.1. DFT, IDFT, and Unnormalized DFT
Discrete Fourier Transform (DFT)


The DFT computes the frequency content of a signal by projecting it onto a basis of harmonic sinusoids:



$$
X = S^H x
$$


In expanded form:



$$
X[k] = \langle x, s_k \rangle = \sum_{n=0}^{N-1} x[n] \cdot \frac{e^{-j\frac{2\pi}{N}kn}}{\sqrt{N}}
$$


    

  • $X[k]$ measures how similar $x[n]$ is to the $k$-th harmonic frequency
    • 

  • This reveals how much of frequency $k$ is present in $x[n]$
    • 





Inverse Discrete Fourier Transform (IDFT)


The IDFT reconstructs the time-domain signal from its frequency components:



$$
x = S X
$$


In expanded form:



$$
x[n] = \sum_{k=0}^{N-1} X[k] \cdot \frac{e^{j\frac{2\pi}{N}kn}}{\sqrt{N}}
$$


This synthesizes $x[n]$ as a linear combination of the harmonic sinusoids weighted by $X[k]$




Unnormalized DFT and IDFT


Some formulations avoid the $\frac{1}{\sqrt{N}}$ normalization by placing the factor entirely in one direction:




    

  • Unnormalized DFT:
    • 


$$
X_u[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j\frac{2\pi}{N}kn}
$$


    

  • Unnormalized IDFT:
    • 


$$
x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{j\frac{2\pi}{N}kn}
$$


This convention is commonly used in numerical implementations (e.g., FFT algorithms)

















3.2. Eigendecomposition and Diagonalization of an LTI System
The Discrete Fourier Transform (DFT) not only decomposes signals into frequency components, but also diagonalizes any linear time-invariant (LTI) system.




Let:


    

  • $H$ be the LTI system matrix (circulant)
    • 

  • $S$ be the DFT matrix (whose columns are harmonic sinusoid eigenvectors)
    • 

  • $\Lambda$ be the diagonal matrix of eigenvalues (frequency response values)
    • 





Then the eigendecomposition of $H$ is:



$$
H = S \Lambda S^H
$$


This expresses $H$ in the frequency domain as a diagonal matrix $\Lambda$, meaning the LTI system simply scales each frequency component independently.




Frequency Response as Eigenvalues


Each diagonal entry of $\Lambda$ corresponds to the system's response at a specific harmonic frequency:



$$
\lambda_k = \sum_{n=0}^{N-1} h[n] \, e^{-j \frac{2\pi}{N} kn} = \langle h, s_k \rangle = H_u[k]
$$


This is the (unnormalized) DFT of the system's impulse response $h[n]$. These values represent the frequency response of the system.




Frequency-Domain Representation


The matrix $\Lambda$ is diagonal, where:



$$
\Lambda =
\begin{bmatrix}
\lambda_0 & & & \\
 & \lambda_1 & & \\
 & & \ddots & \\
 & & &\lambda_{N-1}
\end{bmatrix}
=
\begin{bmatrix}
H_u[0] & 0 & \cdots & 0 \\
0 & H_u[1] & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & H_u[N-1]
\end{bmatrix}
$$




Summary



$$
y = Hx = S \Lambda S^H x
$$
































From right to left:


    

  1. Apply DFT: $S^H x$
    2. 

  2. Multiply by diagonal matrix $\Lambda$
    3. 

  3. Apply inverse DFT: $S (\cdot)$
    4. 

  4. Resulting output: $y$
    5. 





This composition maps the action of a linear time-invariant (LTI) system into the frequency domain, applies a diagonal scaling operation, and transforms the result back to the time domain. It highlights the efficiency and interpretability of the Fourier basis.




In more precise terms:


    

  • The signal $x$ is first transformed into the frequency domain via the DFT matrix $S^H$.
    • 

  • This transformation decomposes $x$ into orthogonal harmonic components, each corresponding to a distinct frequency.
    • 

  • When expressed in this basis, the LTI system becomes diagonal; that is, it scales each frequency component independently without interaction with others.
    • 





This is the essence of diagonalization: a potentially complex linear transformation becomes a simple element-wise multiplication in the frequency domain.


As a result of this decomposition, one avoids analyzing entangled time-domain interactions. Instead, the transformation reveals independent modes of behavior - one per frequency bin.


After applying the desired frequency-domain scaling (via the diagonal matrix $\Lambda$), the inverse DFT (via $S$) reconstructs the output signal in the time domain.




The overall structure is:



$$
y = Hx = S \Lambda S^H x
$$


This relationship offers a powerful interpretation: convolution in the time domain corresponds to multiplication in the frequency domain, and the inverse DFT recovers the time-domain output. This insight simplifies both theoretical analysis and practical implementation of LTI systems.















In [ ]:



    


%% DFT without using the built-in function

x = [1 -1 -1 0 0 0 1 -1];
N = length(x);
X = zeros(1,N);

for k = 0:N-1
    for n = 0:N-1
        X(k+1) = X(k+1) + x(n+1)*exp(-1j*2*pi/N*n*k);
    end
end

X



    















    






X =

  -1.0000 + 0.0000i  -0.4142 + 2.0000i   1.0000 - 0.0000i   2.4142 - 2.0000i   3.0000 + 0.0000i   2.4142 + 2.0000i   1.0000 - 0.0000i  -0.4142 - 2.0000i





















In [ ]:



    


%plot -s 560,600

% plot a discrete-time signal
xn = 0:N-1;
k = 0:N-1;

subplot(3,1,1),
stem(xn,x,'filled'), ylabel('x[n]')
subplot(3,1,2),
stem(k,real(X),'filled'), ylabel('Real of X[k]')
subplot(3,1,3),
stem(k,imag(X),'filled'), ylabel('Imaginary of X[k]')



    















    















    


























In [ ]:



    


%plot -s 560,400

% plot a discrete-time signal

subplot(2,1,1),
stem(xn,x,'filled'), ylabel('x[n]')
subplot(2,1,2),
stem(k,abs(X),'filled'), ylabel('|X[k]|')



    















    















    






























Unnormalized DFT


function [Xk] = dft(xn,N)

% Computes Discrete Fourier Transform
% [Xk] = dft(xn,N)
% Xk = DFT coeff. array over 0 <= k <= N-1
% xn = N-point finite-duration sequence
%  N = Length of DFT
%

n = [0:1:N-1];                       % row vector for n
k = [0:1:N-1];                       % row vecor for k
WN = exp(-1j*2*pi/N);                % Wn factor
nk = n'*k;                           % creates a N by N matrix of nk values
WNnk = WN.^nk;                       % DFT matrix
Xk = xn*WNnk;                        % row vector for DFT coefficients

















$$x[n] = e^{j\frac{2\pi}{8}3n}$$











In [ ]:



    


%plot -s 560,420

k = 3;                      % index for freqeuncy

N = 8;
n = 0:N-1;                  % sampling period

x = exp(1j*2*pi/N*k*n);     % harmonic complex exponential

plot(real(x),imag(x),'o','MarkerFaceColor','b')
grid on
axis equal; axis([-1.5, 1.5,-1.4 1.4])



    















    















    


























In [ ]:



    


%plot -s 560,420

X = dft(x,N);
%X = fft(x,N);

subplot(2,1,1), stem(n,real(x),'filled'), axis tight
xlabel('n'), ylabel('x[n]')
title(['Re(e^{(j2\pi/N kn)}),   k = ',num2str(k)])

subplot(2,1,2); stem(n,abs(X),'filled'), axis tight
xlabel('k')
ylabel('Amplitude of X[k]')
title('Result of FFT')



    















    















    





























$$x[n] = e^{j\frac{2\pi}{8}2n}$$











In [ ]:



    


%plot -s 560,600

% Normalized DFT

k = 2;                        % index for freqeuncy

N = 8;
n = 0:N-1;                    % sampling period

x = exp(1j*2*pi/N*k*n);       % harmonic complex exponential
X = dft(x,N)/N;

subplot(3,1,1), stem(n,real(x),'filled')
xlabel('n'), ylabel('Re\{x[n]\}')
title(['e^{(j2\pi/N kn)},   k = ',num2str(k)],'Fontsize',8)

subplot(3,1,2), stem(n,imag(x),'filled')
xlabel('n'), ylabel('Im\{x[n]\}')
title(['e^{(j2\pi/N kn)},   k = ',num2str(k)],'Fontsize',8)

subplot(3,1,3); stem(n,abs(X),'filled')
xlabel('k','Fontsize',8)
ylabel('Amplitude','Fontsize',8)
title('DFT','fontsize',8)



    















    















    





























$$\text{cos}\;\omega t = \large{e^{i\omega t } + e^{-i\omega t} \over 2}$$





























$$x[n] = \cos \left({\frac{2\pi}{8}1n} \right)$$











In [ ]:



    


%plot -s 560,400

k = 1;                        % index for freqeuncy

N = 8;
n = 0:N-1;                    % sampling period

x = cos(2*pi/N*k*n);          % harmonic complex exponential
X = dft(x,N)/N;

subplot(2,1,1), stem(n,real(x),'filled')
xlabel('n'), ylabel('x[n]'), xlim([0,8])
title(['sin{(2\pi/N kn)},   k = ',num2str(k)],'Fontsize',8)

subplot(2,1,2); stem(n,abs(X),'filled'), ylim([0,1]), xlim([0,8])
xlabel('k','Fontsize',8)
ylabel('Amplitude','Fontsize',8)
title('DFT','fontsize',8)



    















    















    






























In the DFT, the frequency index $k$ maps to the discrete frequency variable $\omega_k = \frac{2\pi}{N}k$. Two common conventions are used for representing this frequency range:




Interval 1: Natural Ordering



$$
0 \leq k \leq N - 1 \quad \Rightarrow \quad 0 \leq \omega_k < 2\pi
$$


This is the standard output ordering of most DFT implementations (e.g., fft in MATLAB or NumPy). Frequencies progress from $0$ up to just below $2\pi$.




Interval 2: Centered (Symmetric) Ordering



$$
-\frac{N}{2} \leq k \leq \frac{N}{2} - 1 \quad \Rightarrow \quad -\pi \leq \omega_k < \pi
$$


This representation centers the zero frequency and is often used for analysis and visualization. It is typically achieved by applying a frequency shift using fftshift.




Note:


The transformation from natural to centered ordering is purely a re-indexing of the same frequency content, but it provides a symmetric view of positive and negative frequencies around $\omega = 0$.















In [ ]:



    


n = 0:N-1;
k = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

X = fft(x,N)/N;
Xs = fftshift(X);

subplot(2,1,1); stem(n,abs(X),'filled'), ylim([0,1]), xlim([0,8])
xlabel('k','Fontsize',8)
ylabel('Amplitude','Fontsize',8)
title('DFT','fontsize',8)

subplot(2,1,2); stem(ks,abs(Xs),'filled'), ylim([0,1]), xlim([-4,4])
xlabel('k','Fontsize',8)
ylabel('Amplitude','Fontsize',8)
title('Shifted DFT','fontsize',8)



    















    















    






























3.3. DFT Properties
The Discrete Fourier Transform (DFT) and its inverse are defined as:



$$
\begin{aligned}
X[k] &= \sum_{n=0}^{N-1} x[n] \cdot \frac{e^{-j \frac{2\pi}{N}kn}}{\sqrt{N}} \\\\
x[n] &= \sum_{k=0}^{N-1} X[k] \cdot \frac{e^{j \frac{2\pi}{N}kn}}{\sqrt{N}}
\end{aligned}
$$


This defines the DFT pair:



$$
x[n] \longleftrightarrow X[k]
$$


DFT Frequencies


    

  • Each $X[k]$ quantifies the similarity between $x[n]$ and the $k$-th harmonic sinusoid:
    • 




$$
s_k[n] = \frac{1}{\sqrt{N}} e^{j \frac{2\pi}{N} kn}
$$


    

  • The index $k$ corresponds to a frequency $\omega_k = \frac{2\pi}{N} k$
    • 





Common Frequency Intervals


Two common conventions are used to interpret the frequency axis of $X[k]$:




    

  • Standard ordering:
    • 




$$
0 \leq k \leq N-1 \quad \Rightarrow \quad 0 \leq \omega_k < 2\pi
$$


This is the default output order of the DFT (e.g., in MATLAB’s fft).




    

  • Centered (symmetric) ordering:
    • 




$$
-\frac{N}{2} \leq k \leq \frac{N}{2} - 1 \quad \Rightarrow \quad -\pi \leq \omega_k < \pi
$$


This form centers the frequency axis around zero and is used for visualization, typically using fftshift.















In [ ]:



    


%plot -s 560,600

%% DFT (FFT) example
N = 16;
h = zeros(1,N);
h(1:8) = 1;
n = 0:N-1;
k = 0:N-1;

X = fft(h,N)/N;
Xs = fftshift(X);

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

subplot(3,1,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(3,1,2),
stem(k,abs(X),'filled'), xlabel('k'), ylabel('abs(X[k])')
subplot(3,1,3),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1])



    















    















    






























DFT and Circular Shift


Let $x[n]$ be a length-$N$ signal, and consider a circular shift of $x[n]$ by $m$ samples:



$$
x[(n - m)_N] \longleftrightarrow e^{-j\frac{2\pi}{N}km} \cdot X[k]
$$


This means:


    

  • The DFT magnitude is unchanged: the amplitude spectrum remains the same.
    • 

  • The DFT phase is modulated by a linear phase term $e^{-j\frac{2\pi}{N}km}$.
    • 





In other words, a circular shift in the time domain corresponds to a linear phase shift in the frequency domain, with no change in spectral content.















In [ ]:



    


%plot -s 800,600

%% DFT and circular shift
N = 16;
h = zeros(1,N);
h(1) = 1;
h2 = cirshftt(h,1,N);

n = 0:N-1;
k = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

X = fft(h,N);
Xs = fftshift(X);

X2 = fft(h2,N);
X2s = fftshift(X2);

subplot(3,2,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(3,2,3),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1])
subplot(3,2,5),
stem(ks,phase(Xs),'filled'), xlabel('k'), ylabel('phase(X[k])'), ylim([-5,5])
xlim([-N/2,N/2-1]), ylim([-5,5])
subplot(3,2,2),
stem(n,h2,'filled'), xlabel('n'), ylabel('x[n]')
subplot(3,2,4),
stem(ks,abs(X2s),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1])
subplot(3,2,6),
stem(ks,phase(X2s),'filled'), xlabel('k'), ylabel('phase(X[k])'), xlim([-N/2,N/2-1]), ylim([-5,5])



    















    















    






























DFT and Modulation


Let $x[n]$ be a length-$N$ signal. Multiplying $x[n]$ by a complex exponential of frequency $\omega_r = \frac{2\pi}{N}r$ results in a frequency shift in the DFT domain:



$$
e^{j\frac{2\pi}{N}rn} \cdot x[n] \longleftrightarrow X[(k - r)_N]
$$


This means:


    

  • Modulation in the time domain (multiplication by a complex exponential) corresponds to a circular shift in the frequency domain.
    • 

  • The DFT is shifted by $r$ frequency bins, with index wrapping via modulo-$N$ arithmetic.
    • 





This duality mirrors the time-shift property:


    

  • Time shift $\rightarrow$ phase shift in frequency
    • 

  • Frequency shift $\rightarrow$ modulation in time
    • 
















In [ ]:



    


%plot -s 800,400

%% DFT and modulation
N = 16;
h = zeros(1,N);
h(1:8) = 1;
h2 = zeros(1,N);

for i = 1:N
    h2(i) = (-1)^i*h(i);
end

n = 0:N-1;
k = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

X = fft(h,N);
Xs = fftshift(X);

X2 = fft(h2,N);
X2s = fftshift(X2);

subplot(2,2,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(2,2,3),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1])
subplot(2,2,2),
stem(n,h2,'filled'), xlabel('n'), ylabel('x[n]')
subplot(2,2,4),
stem(ks,abs(X2s),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1])



    















    















    






























DFT and Circular Convolution




Key property:


Circular convolution in the time domain corresponds to multiplication in the frequency domain.



$$
\begin{align*}
y[n] &= h[n] \otimes x[n] \quad \text{(circular convolution)} \\\\
Y[k] &= H[k]X[k] \quad \text{(pointwise multiplication in frequency)} \\\\
y[n] &= \text{IDFT}(Y[k])
\end{align*}
$$


Proof (Unnormalized DFT):


Let $y[n] = (h \otimes x)[n] = \sum_{m=0}^{N-1} h[(n - m)_N] x[m]$. Then the DFT is:



$$
\begin{align*}
Y_u[k] &=  \sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi}{N}kn} \\\\
&= \sum_{n=0}^{N-1}\left( \sum_{m=0}^{N-1} h[(n-m)_N]x[m]\right) e^{-j\frac{2\pi}{N}kn}\\\\
& = \sum_{m=0}^{N-1}x[m] \left( \sum_{n=0}^{N-1} h[(n-m)_N]e^{-j\frac{2\pi}{N}kn}\right) \\\\
& = \sum_{m=0}^{N-1}x[m] \left( \sum_{r=0}^{N-1} h[r]e^{-j\frac{2\pi}{N}k(r+m)}\right) \\\\
& = \left(\sum_{m=0}^{N-1}x[m]e^{-j\frac{2\pi}{N}km}\right) \left( \sum_{r=0}^{N-1} h[r]e^{-j\frac{2\pi}{N}kr}\right) \\\\
& = X_u[k]H_u[k]
\end{align*}$$


Thus, the unnormalized DFT of the circular convolution is the pointwise product of the unnormalized DFTs:



$$
y[n] = h[n] \otimes x[n] \longleftrightarrow Y_u[k] = H_u[k] X_u[k]
$$




















4. Filtering in Frequency Domain
A fundamental application of the DFT is to implement filtering by manipulating the signal in the frequency domain.




Let $x[n]$ be the input signal, $h[n]$ the impulse response of the system, and $y[n]$ the output signal. Then:



$$
\begin{align*}
\text{Input} \quad &x[n] \\\\
\text{Filter response} \quad &h[n] \\\\
\text{Circular convolution} \quad &y[n] = h[n] \otimes x[n]
\end{align*}
$$


This circular convolution corresponds to multiplication in the frequency domain:



$$
\begin{align*}
X[k] &= \text{DFT of } x[n] \\\\
H[k] &= \text{DFT of } h[n] \\\\
Y[k] &= H[k] \cdot X[k] \\\\
y[n] &= \text{IDFT of } Y[k]
\end{align*}
$$


Compactly:



$$
h[n] \otimes x[n] \quad \longleftrightarrow \quad H[k] \cdot X[k]
$$


The following examples show how filtering operations can be efficiently implemented through element-wise multiplication in the frequency domain.

















Frequency-Domain Interpretation of Filtering


Consider a simple filter with impulse response:



$$
h[n] = [1 \quad 1 \quad 1 \quad 0 \quad \cdots \quad 0]
$$


This filter performs a moving averaging operation over three consecutive samples. It acts as a low-pass filter, meaning it allows low-frequency components of a signal to pass through while attenuating high-frequency components.















In [ ]:



    


%plot -s 560,400

%% low-pass filter
N = 16;
h = zeros(1,N);
h(1:3) = 1;
n = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

X = fft(h,N)/N;
Xs = fftshift(X);

subplot(2,1,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(2,1,2),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2]), ylim([0,0.5])



    















    















    






























In the frequency domain, the DFT of $h[n]$ is concentrated around the low-frequency region. When we multiply $H[k]$ with the spectrum of the input signal $X[k]$, the high-frequency components of $X[k]$ are suppressed, while the low-frequency components are preserved.


As a result, the output signal $y[n]$, obtained by taking the inverse DFT of $Y[k] = H[k] \cdot X[k]$, is a smoothed version of the original input $x[n]$.















In [ ]:



    


%plot -s 560,400

%% high-pass filter
N = 16;
h = zeros(1,N);
h(1) = 1;
h(2) = -1;
n = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);

X = fft(h,N)/N;
Xs = fftshift(X);

subplot(2,1,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(2,1,2),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2]), ylim([0,0.5])



    















    















    


























In [ ]:



    


%plot -s 560,400

%% band-pass filter
N = 16;
h = zeros(1,N);
h(1) = 1;
h(2) = 0;
h(3) = -1;
n = 0:N-1;

kr = [0:N/2-1 -N/2:-1];
ks = fftshift(kr);


X = dft(h,N)/N;
Xs = fftshift(X);

subplot(2,1,1),
stem(n,h,'filled'), xlabel('n'), ylabel('x[n]')
subplot(2,1,2),
stem(ks,abs(Xs),'filled'), xlabel('k'), ylabel('abs(X[k])'), xlim([-N/2,N/2-1]), ylim([0,0.5])



    















    















    


























In [ ]:



    


%plot -s 560,800

N = 50;
n = 0:N-1;
s = sin(pi/N*n) .* [ones(1,N/2), -ones(1,N/2)];

subplot(4,1,1)
stem(s,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal','fontsize',8)

% add noise to the signal
x = s + 0.1*randn(1,N);

subplot(4,1,2)
stem(x,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal + noise','fontsize',8)

% construct moving average filter impulse response of length M
M = 5;
h = ones(1,M)/M;

subplot(4,1,3)
stem(h,'b','Marker','none','LineWidth',2);
axis([1 N 0 1])
title('impulse response','fontsize',8)

% convolve noisy signal with impulse response
y = cconv(x,h,N);

subplot(4,1,4)
stem(y,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('output','fontsize',8)



    















    















    


























In [ ]:



    


%plot -s 800,400

k = n;

H = fft(h,N);
X = fft(x,N);

Y = H.*X;
yi = ifft(Y);

subplot(2,1,1)
stem(n,y,'b','Marker','none','LineWidth',2);
subplot(2,1,2)
stem(n,yi,'b','marker','none','LineWidth',2)



    















    















    


























In [ ]:



    


%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')