7.2.5. Sampling Simulation with MATLAB
The simulation of sampling with MATLAB is complicated by the representation of analog signals and the numerical computation of the analog Fourier transform. Two sampling rates are needed: one being the sampling rate under study, fs, and the other being the one used to simulate the analog signal, fsim >> fs. The computation of the analog Fourier transform of x(t) can be done approximately using the fast Fourier transform (FFT) multiplied by the sampling period. For now, think of the FFT as an algorithm to compute the Fourier transform of a discretized signal.
To illustrate the sampling procedure consider sampling a sinusoid x(t) = cos(2πf0t) where f0 = 1 KHz. To simulate this as an analog signal we choose a sampling period Tsim = 0.5 × 10−4 sec/sample or a sampling frequency fsim = 20,000 samples/sec.
No aliasing sampling—If we sample
x(
t) with a sampling frequency
fs = 6000 > 2
f0 = 2000 Hz, the sampled signal
y(
t) will not display aliasing in its frequency representation, as we are satisfying the Nyquist sampling rate condition.
Figure 7.7(a) displays the signal
x(
t) and its sampled version
y(
t), as well as their approximate Fourier transforms. The magnitude spectrum |
X(Ω)| corresponds to the sinusoid
x(
t), while |
Y(Ω)| is the first period of the spectrum of the sampled signal (recall the spectrum of the sampled signal is periodic of period Ω
s = 2
πfs). In this case, when no aliasing occurs, the first period of the spectrum of
y(
t) coincides with the spectrum of
x(
t) (notice that as a sinusoid, the magnitude spectrum |
X(Ω)| is zero except at the frequency of the sinusoid or ±1 KHz; likewise |
Y(Ω)| is zero except at ±1 KHz and the range of frequencies is [−
fs/2,
fs/2] = [−3, 3] KHz). In
Figure 7.7(b) we show the sinc interpolation of three samples of
y(
t); the solid line is the interpolated values or the sum of sincs centered at the three samples. At the bottom of that figure we show the sinc interpolation, for all the samples, obtained using our function
sincinterp. The sampling is implemented using our function
sampling.
Sampling with aliasing—In
Figure 7.8 we show the case when the sampling frequency is
fs = 800 < 2
fs = 2000, so that in this case we have aliasing. This can be seen in the sampled signal
y(
t) in the top plot of
Figure 7.8(a), which appears as if we were sampling a sinusoid of lower frequency. It can also be seen in the spectra of
x(
t) and
y(
t): |
X(Ω)| is the same as in the previous case, but now |
Y(Ω)|, which is a period of the spectrum of the sampled signal
y(
t), displays a frequency of 200 Hz, lower than that of
x(
t), within the frequency range [−400, 400] Hz or [−
fs/2,
fs/2]. Aliasing has occurred. Finally, the sinc interpolation gives a sinusoid of frequency 0.2 KHz, different from
x(
t).
Similar situations occur when a more complex signal is sampled. If the signal to be sampled is
x(
t) = 2 − cos(
πf0t) − sin(2
πf0t) where
f0 = 500 Hz, if we use a sampling frequency of
fs = 6000 > 2
fmax = 2
f0 = 1000 Hz, there will be no aliasing. On the other hand, if the sampling frequency is
fs = 800 < 2
fmax = 2
f0 = 1000 Hz, frequency aliasing will occur. In the no aliasing sampling, the spectrum |
Y(Ω)| (in a frequency range [−3000, 3000] = [−
fs/2,
fs/2]) corresponding to a period of the Fourier transform of the sampled signal
y(
t) shows the same frequencies as |
X(Ω)|. The reconstructed signal equals the original signal. See
Figure 7.9(a). When we use
fs = 800 Hz, the given signal
x(
t) is undersampled and aliasing occurs. The spectrum |
Y(Ω)| corresponding to a period of the Fourier transform of the undersampled signal
y(
t) does not show the same frequencies as |
X(Ω)|. The reconstructed signal shown in the bottom right plot of
Figure 7.9(b) does not resemble the original signal.
The following function implements the sampling and computes the Fourier transform of the analog signal and of the sampled signal using the fast Fourier transform. It gives the range of frequencies for each of the spectra.
function [y,y1,X,fx,Y,fy] = sampling(x,L,fs)
%
%Sampling
%x analog signal
%L length of simulated x
%fs sampling rate
%y sampled signal
%X,Y magnitude spectra of x,y
%fx,fy frequency ranges for X,Y
%
fsim = 20000; % analog signal sampling frequency
% sampling with rate fsim/fs
delta = fsim/fs;
y1 = zeros(1,L);
y1(1:delta:L) = x(1:delta:L); y = x(1:delta:L);
% analog FT and DTFT of signals
dtx = 1/fsim;
X = fftshift(abs(fft(x))) ∗ dtx;
N = length(X); k = 0:(N− 1); fx = 1/N.*k; fx = fx ∗ fsim/1000− fsim/2000;
dty = 1/fs;
Y = fftshift(abs(fft(y))) ∗ dty;
N = length(Y); k = 0:(N− 1); fy = 1/N.*k; fy = fy∗ fs/1000− fs/2000;
The following function computes the sinc interpolation of the samples.
function [t,xx,xr] = sincinterp(x,Ts)
%
% Sinc interpolation
% x sampled signal
% Ts sampling period of x
% xx,xr original samples and reconstructed in range t
%
N = length(x)
t = 0:dT:N;
xr = zeros(1,N ∗ 100 + 1);
for k = 1:N,
xr = xr + x(k)∗ sinc(t−(k− 1)); end
xx(1:100:N ∗ 100) = x(1:N);
xx = [xx zeros(1,99)];
NN = length(xx)
t = 0:NN− 1;t = t ∗ Ts/100;