Discrete Signal Processing (DSP):

Discrete Fourier Transformation (DFT)

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Eigenanalysis

1.1. Eigenvector and eigenvalues

• Given matrix $A$,

$$ Av=\lambda v$$

• Using this, we can find the following equation

$$ \begin{align*} AV &=\, V\Lambda\\\\ AV &= \, [v_0 \mid v_1 \mid \cdots \mid v_{N-1}] \begin{bmatrix} \lambda_0 & & & \\ & \lambda_1 & & \\ & & \ddots & \\ & & &\lambda_{N-1} \end{bmatrix} \end{align*}$$

• We can change to

$$ \begin{align*} A &= V\Lambda V^{-1} \implies \text{Eigen-decomposition}\\ \\ V^{-1} AV\, &= \, \Lambda \end{align*} $$

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

• For length-$N$ signals, $H$ is an $N \times N$ circulent matrix with entries

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

$\quad\;$ where $h$ is the impulse response

• Goal: Calcuate the eigenvectors and eigenvalues of $H$

• Eigenvectors $v$ are input signals that emerge at the system output unchanged (except for a scaling by the eigenvalue $\lambda$) and so are somehow "fundamental" to the system

• Fact: the eigenvectors of a circulent matrix $H$ for a LTI system are the complex harmonic sinusoids

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

• Prove that harmonic sinusoids are the eigenvectors of LTI systems simply by computing the circular convolution with input $s_k$ and applying the periodicity of the harmonic sinusoids.

$$\begin{align*} 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*}$$

• The eigenvalue $\lambda \in \mathbb{C}$ corresponding to the sinusoid eigenvector $s_k$ is called the frequency response at frequency $k$ since it measures how the system "responds" to $s_k$

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

• $\lambda _k$ means the number of $s_k$ in $h[n]$ $\implies$ similarity

1.3. Eigenvector Matrix of Harmonic Sinusoids

Stack $N$ normalized harmonic sinusoid $\{s_k\}_{k=0}^{N-1}$ as columns into an $N\times N$ complex orthonormal basis matrix

$$ \begin{align*} S &= \left[s_0 \mid s_1 \mid \cdots \mid s_{N-1} \right] \\ \\ s_k[n] &= \frac{1}{\sqrt{N}}e^{j\frac{2\pi}{N}kn} \end{align*} $$

%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 by Harmonic Sinusoids

• We now turn back to linear algebra to understand transforms, which map signals between different "domains"

• Recall that signals can be interpreted as vectors in a vector space

• We now review the concepts of a basis for a vector space

• As we will see, different signal transforms (and "domains") correspond to different basis

2.1. Basis

A basis $\{b_k\}$ for a vector space $V$ is a collection of vectors from $V$ that linearly independent and span $V$

• Span: all vectors in $V$ can be represented as a linear combination of the basis vectors $\{b_k\}_{k=0}^{N-1}$

$$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}$$

• Linearly independent: none of the basis vectors can be represented as a linear combination of the other basis vectors

• Basis matrix: stack the basis vectors $b_k$ as columns

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

• Using this matrix $B$, we can now write a linear combination of basis elements as the matrix/vector product

$$ \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*}$$

2.2. Orthonormal Basis

• An orthogonal basis $\{b_k\}_{k=0}^{N-1}$ for a vector space $V$ is a basis whose elements are mutually orthogonal

$$\langle b_k,b_l\rangle \, =\, 0\, \quad k\, \neq\, l$$

• An orthonormal basis$\{b_k\}_{k=0}^{N-1}$ for a vector space $V$ is a basis whose elements are mutually orthogonal and normalized (in the 2-norm)

$$\begin{align*} \langle b_k,b_l\rangle \, &=\, 0\, \quad k\, \neq\, l, \quad \text{ and}\\\\ \| b_k\|_2 \, &= \, 1 \quad \forall k \end{align*}$$

• Then $B$ is a unitrary matrix

$$B^H B = I \implies B^{-1} = B^H, \qquad \text{where } B^H \text{ is Hermitian (complex conjugate) transpose}$$

$$ \begin{align*} B^H B = \begin{bmatrix} b_0^H\\ b_1^H\\ \vdots\\ b_{N-1}^H \end{bmatrix} \left[b_0 \mid b_1 \mid b_2 \mid \cdots \mid b_{N-1} \right] \end{align*} = \begin{bmatrix} 1 & \\ & 1\\ & & \ddots\\ & & & 1 \end{bmatrix} $$

2.3. Signal Representation by Orthonormal Basis

Given an orthonormal basis $\{b_k\}_{k=0}^{N-1}$ and orthonormal basis matrix $B$, signal representation for any signal $x$

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

• Synthesis: Build up the signal $x$ as a linear combination of the basis elements $b_k$ weighted by the weights $\alpha_k$

• Analysis: Compute the weights $\alpha_k$ such that the synthesis produces $x$; the weight $\alpha_k$ measures the similarity between $x$ and the basis element $b_k$

2.4. Harmonic Sinusoids are an Orthonormal Basis

Stack $N$ normalized harmonic sinusoid $\{s_k\}_{k=0}^{N-1}$ as columns into an $N\times N$ complex orthonormal basis matrix

$$ \begin{align*} S = \left[s_0 \mid s_1 \mid \cdots \mid s_{N-1} \right] \quad \text{where}\quad s_k[n] = \frac{1}{\sqrt{N}}e^{j\frac{2\pi}{N}kn} \end{align*} $$

$$S^H S = I \implies \langle s_k,s_l\rangle \, =\, 0 \quad k\, \neq l,\quad \text{and}\quad \lVert s_k \rVert_2\, =\, 1$$

3. Discrete Fourier Tranformation (DFT)

Jean Baptiste Joseph Fourier had the radical idea of proposing that all signals could be represented as a linear combination of sinusoids.

3.1. DFT, IDFT, Unnormalized DFT

DFT

• The weight $X[k]$ measures the similarity between $x$ and the harmonic sinusoid $s_k$

• It finds the “frequency content” of $x$ at frequency $k$

$$ \begin{align*} X\, &=\, S^Hx \\\\ X[k]\, &=\, \langle x,s_k\rangle\, =\, \sum \limits_{n=0}^{N-1} x[n]\,\frac{e^{-j\frac{2 \pi}{N}kn}}{\sqrt{N}} \end{align*}$$

Inverse DFT

• It is returning to time domain
• It builds up the signal $x$ as a linear combination of $s_k$ weighted by the $X[k]$

$$ \begin{align*} x\, &=\, SX \\\\ x[n]\, &=\, \sum \limits_{k=0}^{N-1}X[k]\, \frac{e^{j\frac{2 \pi}{N}kn}}{\sqrt{N}} \end{align*} $$

Unnormalized DFT

• Normalized forward and inverse DFT

$$ \begin{align*}X[k] &= \sum \limits_{n=0}^{N-1} x[n]\,\frac{e^{-j\frac{2 \pi}{N}kn}}{\sqrt{N}}\\\\ x[n] &= \sum \limits_{k=0}^{N-1}X[k]\, \frac{e^{j\frac{2 \pi}{N}kn}}{\sqrt{N}} \end{align*} $$

• Unnormalized forward and inverse DFT

$$ \begin{align*}X_u[k] &= \sum \limits_{n=0}^{N-1} x[n]\,e^{-j\frac{2 \pi}{N}kn}\\\\ x[n] &= \frac{1}{N}\sum \limits_{k=0}^{N-1}X[k]\, e^{j\frac{2 \pi}{N}kn} \end{align*} $$

3.2. Eigendecomposition and Diagonalization of an LTI Systems

• $H$ is LTI System matrix

• $S$ is harmonic sinusoid eigenvectors matrix

• $\Lambda$ is eigenvalue diagonal matrix

• The eigenvalues are the frequency response (unnormalized DFT of the impulse response)

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

• Place the $N$ eigenvalues $\{\lambda_k\}_{k=0}^{N-1}$ on the diagonal of an $N \times N$ matrix

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

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

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

%% 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

%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]')

%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}$$

%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])

%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}$$

%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)$$

%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)

Typical interval 1: $0\leq k \leq N-1$ corresponds to frequencies $\omega_k$ in the interval $0 \leq \omega \leq 2\pi$

Typical interval 2: $-\frac{N}{2}\leq k \leq \frac{N}{2}-1$ corresponds to frequencies $\omega_k$ in the interval $-\pi \leq \omega \leq \pi$

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

$$ \begin{align*} X[k]\, &= \sum \limits_{n=0}^{N-1} x[n]\,\frac{e^{-j\frac{2 \pi}{N}kn}}{\sqrt{N}}\\ x[n]\, &=\, \sum \limits_{k=0}^{N-1}X[k]\, \frac{e^{j\frac{2 \pi}{N}kn}}{\sqrt{N}} \end{align*}$$

DFT pair

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

DFT Frequencies

• $X[k]$ measures the similarity between the time signal $x[n]$ and the harmonic sinusoid $s_k[n]$

• $X[k]$ measures the "frequency content" of $x[n]$ at frequency $\omega_k = \frac{2\pi}{N} k$

• Typical interval 1: $0 \leq k \leq N-1$ corresponds to frequencies $\omega_k$ in the interval $0 \leq \omega \leq 2\pi$

• Typical interval 2: $-\frac{N}{2} \leq k \leq \frac{N}{2} -1$ corresponds to frequencies $\omega_k$ in the interval $-\pi \leq \omega \leq \pi$

%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

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

• no amplitude changed

• phase changed

%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

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

%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

• Circular convolution in the time domain $=$ multiplication in the frequency domain

$$ \begin{align*} Y[k] &= H[k]X[k] \\\\ h[n]\otimes x[n] &\longleftrightarrow H[k]X[k] \\\\ y[n] &= \text{IDFT}(Y[k]) \end{align*} $$

• Proof

$$ \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*}$$

4. Filtering in Frequency Domain

$$ \begin{align*} Y[k] &= H[k]X[k] \\\\ h[n]\otimes x[n] &\longleftrightarrow H[k]X[k] \\\\ y[n] &= \text{IDFT}(Y[k]) \end{align*} $$

%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])

%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])

%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])

%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)

%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)

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