% EE-3220 Digital Signal Processing
% Monday 15 December 2014
% Dr. Durant
function l02alias
% Try changing the values of f (frequencies of sinusoidal signal) below. You can have 1 or more, but they must be in a column vector.
% Also try changing the value of fs (sampling frequency) below to demonstrate aliasing.
f = [400; 700]; % Hz, column vector, cosine frequencies
tmax = 0.005; % End time for graphs, 2/min(f) is reasonable (gives 2 cycles of slowest wave)
t = linspace(0,tmax); % time from 0 to tmax, default 100 samples for plenty of detail
yanlg = cos(2*pi*f*t); % cosine for each f at all times t (1 row for each f, note that f*t is an OUTER product)
figure, plot(t,yanlg,'b-') % Plot cosines vs. time, blue lines
xlabel('Time (s)'), ylabel('Signal voltage (V)')
fs = 1000; % sampling frequency, Hz
N = ceil(tmax*fs); % number of samples
tSamp = linspace(0, (N-1)/fs, N); % sampling time: start, finish, number of samples
hold on % add to existing plot instead of replacing it
for tv = tSamp % foreach sampling time,
% plot circles representing sampled value for all frequencies...
plot(tv, cos(2*pi*f*tv), 'ko'); % black circles
end
% For each frequency component, report whether it is aliased, and draw the aliased continuous time waveform for the aliased components.
omega = f/fs * 2*pi; % convert analog to digital frequency
for idx=1:length(f) % for each frequency
if abs(omega(idx)) <= pi
fprintf('|omega| = %g * pi <= pi, so no aliasing occurs for the %g Hz component.\n', ...
abs(omega(idx))/pi, f(idx))
else
omegaAlias = mod(omega(idx),2*pi); % adjust by 2k*pi to make range [0,2pi)
if omegaAlias > pi % adjust range to [-pi,pi)
omegaAlias = omegaAlias - 2*pi; % e.g., 3/2 pi becomes -1/2 pi, always want the minimum magnitude digital frequency ("baseband assumption" throughout EE3220)
fAlias = omegaAlias/(2*pi)*fs; % convert to Hz
end
fprintf('|omega| = %g * pi > pi, so aliasing to %g * pi occurs.\nAt fs = %g, that is %g Hz instead of the %g Hz in the original signal.\n', ...
abs(omega(idx))/pi, omegaAlias/pi, fs, fAlias, f(idx))
plot(t,cos(2*pi*fAlias*t),'r-') % plot the aliased analog waveform; this will come out of the reconstruction filter (system/reconstruction delay not included here)
end
end
end % function