Assignment title: Information
1
EEE 304 Lab Exercise 2: Sampling and Reconstruction
1. Sampling Theory
Suppose we sample a signal, such as the one in Fig. 1, which is band-limited to frequencies between 0 Hz
and fn Hz. If we want to reconstruct the signal later, we must ensure that the sample rate, fs, is strictly
greater than 2*fn. A signal sampled at fs = 2*fn is said to be Nyquist sampled, and fs = 2*fn is called the
Nyquist frequency. This is a consequence of the Nyquist-Shannon Sampling Theorem which states that in
order to perfectly capture or reconstruct a signal after sampling it, the sampling rate must be at least
twice that of the highest frequency (fn) in the signal band.
Note that information may be lost if a signal is sampled exactly at the Nyquist frequency. For example, the
sine wave in Fig. 1 has a frequency of 1 Hz. The Nyquist frequency is therefore 2 Hz. If we sample the sine
wave at a rate of 2 Hz, say at t=0, t=0.5, t=1, and so on, all the sample values will be 0. The signal will look
as if it were identically 0, and no reconstruction method will be able to recreate the 1 Hz sine wave. This is
why fs must be strictly greater than 2*fn.
Fig 1. Sine Wave Signal
2. Reconstruction Theory
The reconstruction of a signal is implemented by applying a low pass filter to the sampled signal. As an
example, consider person 'A' speaking into a telephone. Typical human voice signals are between 300 Hz -
3.4 kHz. A sampling rate of 8 kHz is more than sufficient to accurately capture all the signals present in
person 'A's voice. The signal is digitized and sent over a communication protocol to person 'B'. A low pass
filter filters the 8 kHz signal back to 0.3 - 3.4 kHz and accurate reconstruction of person 'A's voice is heard
through the speaker by person 'B'.
3. Frequency-Magnitude Character of Signal
The frequency-magnitude (or spectrum) of discrete signal is described by the Discrete Fourier Transform
(DFT). In MATLAB, the DFT is done be the command 'fft'. Use the 'help fft' in Matlab to learn how to use it.
You may also refer to Lab 1 Manual.
2
The relationship between the normalized analog frequency, ω, and digital frequency, f, is defined below
ω = 2πf/fs
where, fs is the sampling frequency. Therefore, ω =2π corresponds to f=fs. For 100 points fft, these 100
points occupy the range, [0, 2π] evenly. The frequency interval between 2 consecutive points is fs/(100-1).
The sample code below shows a way to plot the double-sided frequency-magnitude plot.
%% let's start with a clean slate
clear all % clears all workspace variables
close all % closes all windows and figures
clc % clears the command window
%% create the sine wave
fs=1000; % sample frequency 1000Hz
t=0:1/fs:0.2; % time duration: 0.2 second with sampling
period 1/fs
fn=50; % single frequency 50Hz
y=sin(2*pi*fn*t); % the sine wave
figure % create a figure
plot(t,y); % plot the wave
title('the sine wave'); % title of the figure
xlabel('t (second)'); % the label of x-axis
ylabel('y(t)'); % the label of y-axis
%% plot the frequency response. Double sided spectrum
no_pts = 2^nextpow2(length(y)); % create the no. of points of FFT
Y = fft(y,no_pts)/(no_pts); % generate FFT and divide by no_pts to
normalize fft ;
f = fs*linspace(0,1, length(Y) ) - fs/2 ; % shift to be double-sided spectrum
figure
plot( f(length(Y)/2+1:end),abs(Y(1:length(Y)/2)), 'b' ), hold on %plot positive
frequency
plot( f(1:length(Y)/2+1),abs(Y(length(Y)/2:end)), 'b' ); %plot negative
frequency
ylabel('|Y(f)|'); xlabel('f (Hz)')
figure
plot( f, fftshift(abs(Y)) ) % fftshift generates the same result as the two
plot commands shown above
4. Some Derivations: Continuous Fourier Transform (CFT) interpretation of the
Discrete Fourier Transform (DFT), the Fast Fourier Transform (FFT), and the
Discrete Time Fourier Transform (DTFT)
The DTFT is a discrete-time transform defined for an infinite sequence of numbers, as follows.
n
j j n
X DTFT (e ) x(n)e
The DTFT is periodic with period 2π.
The DFT is a discrete-time transform defined for a finite sequence of N numbers, as follows.
3
1 0
2
( ) ( )
N n
kn
N
j
X DFT k x n e
The term FFT is used to signify a specific method to compute the DFT in an efficient manner. While a
transform in its own right, the DFT/FFT can be interpreted in terms of the CT-Fourier transform (CTF) for a
sampled signal x(t). Let us consider a signal x(t), sampled at the time instants nT, n = 0,1,…N-1.
1) Periodic Case:
We assume that the signal is periodic both in discrete and continuous time and that the period is NT.
The CFT is found by computing the Fourier Series (FS) expansion first:
NT
x t e dt w
NT
X jw a w kw a
NT
jkw t
k
k
CFT k
2
( ) 2 ( 0) ; 1 ( ) 0 ; 0
The DTFT is also found by computing the DFS expansion first:
N
x n e
N
X e b k b
N
jk n
k
k
k
j
DTFT
2
( ) 2 ( 0) ; 1 ( ) 0 ; 0
Notice
k
N
jkw nT
k x nT e T a
NT
x n x nT b
( ) ( ); 1 ( ) 0
The last approximation is the discretization of the FS integral.
The DFT is:
k
N
jk n
X DFT k x n e Nb
( ) ( ) 0 (1)
If, in addition, the signal is bandlimited and its maximum frequency is smaller than half of the sampling
frequency, π/T, then there is no aliasing and ak bk .
2) Aperiodic Case:
Suppose that the signal is 0 outside the sampled interval [0, (N-1)T]. The Fourier transform of this signal is
N T
jwt jwt
X CFT jw x t e dt x t e dt
( 1)
0
( ) ( ) ( )
First, in a naïve approach, let us approximate the integral by an Euler discretization, whereby the
integrand is taken as piecewise constant.
1 0
( 1)
0
( ) ( ) ( )
N n
jwnT
N T
jwt
X CFT jw x t e dt x nT e T
Let us consider now a sampled version of the Fourier transform at the frequencies
TN
w kw k 2
0 ,
where k=0,1,…,N-1. Then,
1 0
1 2
0
2
( 0) ( ) ( )
N n
kn
N
N j
n
nT
NT
k
j
X CFT jkw x nT e T T x nT e
Thus,
CFT x t DFT x nT T FFT x nT T
{ ( )}|w kw { ( )} { ( )}
0
(2)
This approximation is valid as long as the integrand approximation by a piecewise constant function is
reasonable. This implies that T must be small and x(t) should not change significantly inside any sampling
interval of length T.
4
In a more precise formulation, the sampled signal
k
x(t) (t nT ) has a Fourier transform
k
s
k
CFT s s X j w kw
T
w kw
T
X jw X ( jw)* 2 ( ) 1 ( ( ))
1 2
( )
_
where, as usual,
T
w
s
2
. This is the familiar form of the sampled signal Fourier transform as a
summation of shifted replicas of the signal Fourier transform. Of course, aliasing effects will occur if XCFT(jw)
extends beyond half the sampling frequency. On the other hand, using the Fourier transform definition on
the sampled signal we find
n
jwnT
n
jwnT
n
jwnT
jwt
k n
CFT s
x nT e
x nT t nT e dt x nT e t nT dt
X jw F x t t nT x t t nT e dt
( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( )
_
Again, evaluating this transform at a discrete set of frequencies
TN
w kw k 2
0 , we get
X jkw x nT e DFT x nT FFT x nT CFT x t w kw T
n
kn
N
j
CFT s ( ) { ( )} { ( )} { ( )}| /
0
2
_ 0
(3)
In other words, the FFT of the sampled signal is equal to the Fourier transform of the sampled signal
evaluated at the frequencies kw0, in the interval [0, ws]. It will also be approximately equal to XCFT(jw)/T (at
the corresponding discrete frequencies in the interval [0, ws/2]) as long as any aliasing effects are small.
Finally, in the aperiodic case, a comparison of the DFT and DTFT definitions
n
j j n
X DTFT (e ) x(n)e
1 0
2
( ) ( )
N n
kn
N
j
X DFT k x n e
yields that as long as x(n) is zero outside [0,N-1],
( ) ( )
2 /
X e X k
DFT
k N
j
DTFT
(4)
Note there will be some aliasing due to the fact that a finite duration signal cannot be bandlimited. This
implies that XCFT(jw) is only approximately equal to the replicas in XCFT_s(jw).
5. Derivation of Fourier Transform of a Windowed Sinusoid
In this lab exercise, we define the signal y(t) = sin(w0t), for t in [0,Tmax]. Consequently, the CFT of y(t) is
2
( ) /
0
( ) / 2 0 max
0
0 max
max / 2
0 [0, ] 0 0
0 max 0 max
max
max
1 sin(( ) / 2) sin(( ) / 2)
2sin( / 2)
( ) ( ) *
1 2
{sin( ) ( )}
j w w T j w w T
jwT
T
e
w w
w w T
e
w w
w w T
j
e
w
wT
w w w w
j
F w t pulse t
5
(5)
This expression is used as a basis to compare the FFT-computed transform. The latter will include aliasing
effects, attaining its theoretical value as the sampling interval approaches zero.
The following MATLAB code compares the CFT of a windowed sinusoid derived by Equation (5) and the
FFT method:
close all; clear all;
%%
% plot the time domain signal
% Define the sampling frequency and the time vector. (8000 points / 44 kHz ~
0.2 s worth of data)
% Plot the output and look at the first 0.004s. Play the sound. you must
% have a speaker or headphone connected to listen
fs=44100; no_pts=8192;
t=([0:no_pts-1]')/fs;
y1=sin(2*pi*1000*t);
figure;
plot(t,y1);
xlabel('t (second)')
ylabel('y(t)')
axis([0,.004,-1.2,1.2]) % constrain axis so you can actually see the wave
sound(y1,fs); % play sound using windows driver.
%%
% Check the frequency domain signal. fr is the frequency vector and f1 is the
magnitude of F{y1}.
fr=([0:no_pts-1]')/no_pts*fs; %in Hz
fr=fr(1:no_pts/2); % single-sided spectrum
f1=abs(fft(y1)); % compute fft
f1=f1(1:no_pts/2)/fs;
%%
% F is the continuous time Fourier. (See derivation notes.) Notice the small
% amount of aliasing due to the fact that the truncated sinusoid is not
bandlimited.
frp=fr*2*pi;tmax=max(t);
F1=1/1i*sin((frp-1000*2*pi)*tmax/2).*exp(-1i*(frp-1000*2*pi)*tmax/2)./(frp-
1000*2*pi); % li = i, the imaginary part symbol
F2=1/1i*sin((frp+1000*2*pi)*tmax/2).*exp(-
1i*(frp+1000*2*pi)*tmax/2)./(frp+1000*2*pi);
F=abs(F1-F2); % magnitude
figure;
plot(fr, F, fr, f1) % compare the continuous time Fourier with FFT, linear
scale
axis([0,1500,0,0.09]) % constrain axis so you can actually see the pulse
xlabel('f (Hz)')
ylabel('|Y(f)|')
legend('Continuous-time FT', 'FFT')
figure;
loglog(fr,F,fr,f1); % compare the continuous time Fourier with FFT, log-log
scale
xlabel('f (Hz)')
ylabel('|Y(f)|')
legend('Continuous-time FT', 'FFT')
6. Aliasing
6
Now we sample the same sinusoid (1 kHz) in Section 5 at a lower sampling rate to observe the aliasing
effects. You must have all the variables from the code in Section 5 in the MATLAB workspace for this to
work. Also observe the impact of zero-order hold (ZOH) signal reconstruction on the quality of the output
signal.
close all;
%%
% Sample the sinusoid at 44/4 kHz. Compare the two output and play the sound.
% This may still sound OK because of the internal filtering of the soundcard.
a=4; t_a=([0:no_pts/a-1]')/fs*a;
y_a=sin(2*pi*1000*t_a);
figure
plot(t,y1,t_a,y_a);
xlabel('t (second)')
ylabel('y(t)')
legend('fs = 44kHz', 'fs = 11kHz')
axis([0,0.004,-1.2,1.2]);
display('sampled at 11kHz,played at 11kHz')
sound(y_a,fs/a); % specify sampling rate for the soundcard
pause(2) % pause for 2 sec to separate from the next sound command
%%
% Check the frequency domain signal. Notice that the replicas of F{y_a} are
% now in the audible range (not shown).
fr_a=([0:no_pts/a-1]')/no_pts*a*fs/a;
fr_a=fr_a(1:no_pts/a/2);
f_a=abs(fft(y_a));
f_a=f_a(1:no_pts/a/2)/fs*a;
figure
loglog(fr,F,fr,f1,fr_a,f_a); % log-log scale
xlabel('f (Hz)')
ylabel('|Y(f)|')
legend('Continuous-time FT', 'FFT (fs=44kHz)', 'FFT (fs=11kHz)')
%%
% Signal reconstruction:
% Use the interp1 function to get the ZOH (zero-order hold) version of the signal at 44 kHz.
% The ZOH version is how the signal would sound with a 44 kHz D/A converter
% soundcard
% Its difference from the original is significant in both the time and frequency domain and its reproduction is quite
poor.
y_n = interp1(t_a,y_a,t,'nearest','extrap');
f_n=abs(fft(y_n));
f_n=f_n(1:no_pts/2)/fs;
figure
plot(t,y1,t,y_n);
axis([0,0.004,-1.2,1.2]);
xlabel('t (second)')
ylabel('y(t)')
legend('fs = 44kHz', 'fs = 11kHz reconstructed at 44kHz using ZOH')
figure
loglog(fr,F,fr,f1,fr,f_n);
xlabel('f (Hz)')
ylabel('|y(f)|')
legend('Continuous-time FT', 'FFT (fs = 44kHz)', 'FFT (fs = 11kHz reconstructed at 44kHz using ZOH)')
display('sampled at 44kHz,played at 44kHz')
sound(y1,fs);
pause(2)
display('sampled at 11kHz,played at 44kHz')
7
sound(y_n,fs);
7. Upsampling
Given a discrete-time signal y(n), the upsampled signal (by a) is defined as the signal z(n) such that
z(n) y( 0k) otherwise if n ak
Equivalently, we may think of z(n) as the signal produced by sampling the continuous time signal y(t) with
the impulse train 𝑝𝑢(𝑡) = ∑ 𝑏𝑛 𝛿(𝑡 − 𝑛𝑇
𝑎
), where bn = 1, if n/a is an integer, and 0 otherwise.
Clearly, the upsampling impulse train is identical with the original delta train 𝑝(𝑡) = ∑ 𝛿(𝑡 − 𝑛𝑇), and
therefore, the continuous-time sampled signals y(t)p(t) and y(t)pu(t) are the same, and have the same
Fourier transforms. The upsampled signal, however, has a smaller sampling time (larger sampling
frequency), implying that the FFT of the discrete-time sequence will contain several replicas of F{y}. The
upsampling process does not create any new information about the time signal but allows for the
implementation of better low-pass filters in discrete time.
Instead of using ZOH to reconstruct the signal as shown in Section 6, the following MATLAB code use
upsampling to "interpolate" the values of the low rate signal and convert it to a high rate signal, and then
pass the signal through a low pass filter to recover the original signal. This does not overcome any aliasing
effects that occurred at sampling but improves the reproduction by reducing the undesirable properties
of the low rate ZOH. Note this implementation does require a better high rate D/A converter. In the
MATLAB code below, a digital filter is implemented by approximating the impulse response of an ideal
low-pass filter. Once again, variables from the previous experiments must be present in MATLAB
workspace for this script to run correctly.
close all;
%%
% Generate an "Upsampled" version of the low-rate signal at 44 kHz.
% zero-fill the added sampling points
% The upsampled signal contains several replicas of the original Fourier
% transform within the sampling frequency and it is not an approximation
% of the original signal.
y_u=y1*0;k=1:a:no_pts;y_u(k)=y_a;
f_u=abs(fft(y_u));
f_u=f_u(1:no_pts/2)/fs*a;
figure
plot(t,y1,t,y_u);axis([0,0.004,-1.2,1.2])
xlabel('t (second)')
legend('y(t)','y(t) upsampled')
figure
plot(fr,f1,fr,f_u, '--'); % linear scale
xlabel('f (Hz)')
ylabel('|Y(f)|')
legend('FFT (fs = 44kHz)', 'FFT (fs = 11kHz, upsampled)')
%%
% However, the original signal can be "recovered" after low-pass filtering.
% Here, the filter is defined in terms of its impulse response and the output
% is computed as a convolution sum. Notice that the impulse response HS is 1999
% points long yielding a 999 points delay in the sound reproduction. This is
% significant but not necessarily impractical.
% The low pass filter is a sinc function, whose frequnce response is an
% ideal lowpass filter with cutoff frequency of 3kHz.
hs=sin(3000*2*pi*t)/pi./t; i=1000:-1:2;
%make it symmetric, sinc(0)=sin(2*pi*f*t)/(2*pi*f*t)=1, so
%sin(2*pi*f*0)/(pi*0)=2*f
HS=[hs(i);3000*2;hs(2:1000)]*a/fs;
8
figure
plot(HS); % see the response of the filter designed
xlabel('t')
ylabel('h(t)')
y_f=conv(HS,y_u); % compute the filtered signal
%figure;
%plot(y_f) % uncomment to see the filtered signal
y_f=y_f(1000:999+no_pts); % extract the middle part of the signal
figure
plot(t,y1,t,y_f);axis([0,0.004,-1.2,1.2])
xlabel('t (second)')
legend('y(t)','y(t) upsampled-filtered')
f_f=abs(fft(y_f));f_f=f_f(1:no_pts/2)/fs; % get fft of y_u
figure
loglog(fr,f1,fr,f_f); % log-log scale
xlabel('f (Hz)')
ylabel('|Y(f)|')
legend('FFT (fs = 44kHz)', 'FFT (fs = 11kHz, upsampled-filtered)')
disp('sampled at 44kHz, played at 44kHz');
sound(y1,fs);
pause(2)
disp('sampled at 11kHz, upsampled-filtered, played at 44kHz');
sound(y_f,fs);
8. Equalizing
Equalizers are often used in audio equipment to amplify (or attenuate) frequency bands in an effort to
improve the quality of reproduction of the original signal. Their objective is to cancel (or invert) the audio
signal distortion caused by equipment limitations (amplifier, speaker, media) or the environment where
the signal is reproduced.
The following MATLAB code uses digital low-pass filters to implement a "frequency equalizer." The test
signal is the standard Windows sound "Tada." The MATLAB program loads the sound (make sure its copy
is in the MATLAB working directory) and filters it with three different low-pass filters. Then, by taking
different linear combinations of the resulting signals we can amplify or attenuate the energy of the signal
in the desired frequency band(s).
When you run the program, it will ask for equalizer coefficients that it uses to scale up the signal
components in four frequency bands. Enter a vector such as [1 100 100 1] will scale the signal in the
middle two frequency bands by 100, and it is essentially a bandpass filter.
% Equalization
clear all;close all;clc;
[y,f]= audioread('tada.wav');
y=y(:,1); % keep only one channel
% plot fft of the sound signal
ly=length(y); % no of points in the data
ty=([1:ly]'-1)/f;
w=([1:ly]'-1)/ly*f;w=w(1:ly/2); %in Hz, single-sided spectrum
fy=abs(fft(y));fy=fy(1:ly/2)/f; % FFT of original signal
% define 3 999-point lowpass filters
i=2:500;i_=500:-1:2; L=500:499+ly;
hs1=sin(700*2*pi*ty)/pi./ty;HS1=[hs1(i_);700*2;hs1(i)]/f; % cutoff freq = 700Hz
hs2=sin(2000*2*pi*ty)/pi./ty;HS2=[hs2(i_);2000*2;hs2(i)]/f; % cutoff freq =
2kHz
hs3=sin(5000*2*pi*ty)/pi./ty;HS3=[hs3(i_);5000*2;hs3(i)]/f; % cutoff freq =
5kHz
z1=conv(HS1,y);z1=z1(L);
9
z2=conv(HS2,y);z2=z2(L);
z3=conv(HS3,y);z3=z3(L);
done=0;
coef=input('enter equalization coefficients [low, mid-low, mid-hi, hi]: ');
if isempty(coef),coef=[1 1 1 1];end
% extract the signal in each of the frequency band: 0~700, 700~2k,
% 2k~5k, >5k, multiply each with the corresponding equalizer coeffs
% entered and add up the 3 scaled signal
yx=[z1 z2-z1 z3-z2 y-z3]*coef';
fx=abs(fft(yx));fx=fx(1:ly/2)/f;
figure
loglog(w, fy, w,fx, '--');
legend('orignal', 'equalized signal')
xlabel('f (Hz)');
ylabel('|Y(f)|');
disp('original')
soundsc(y,f)
pause(2)
disp('equalized')
soundsc(yx,f)