% Simple script to illustrate continuous-time sampling with aliasing
% Last revised 9/7/99 - rms



pause
axis on
fs = 10000;	% sampling rate in Hz
t = 0:(1/fs):.5; %10 kHz sampling rate
fhi = 44100;	% CD sampling rate
thi  = 0:(1/(44100)):.5; 
tfast =  0:(1/(10*fs)):.001;	% Only show first 1 msec
f0 = 1000;	%1000 Hz
w0 = 2*pi*f0; 
fnyq = fs/2;


for harmonic  = 1:9,
x = cos(harmonic*w0*t);
xhi = cos(harmonic*w0*thi);
harmsamp = min(harmonic,((2*fnyq/f0)-harmonic));

xfast = cos(harmonic*w0*tfast);
yfast = cos(harmsamp*w0*tfast);

subplot(2,1,1),plot(tfast,xfast,'r')
hold on
stem(t(1:11),x(1:11))
hold off
str1 = num2str(harmonic);
title(['Input frequency ',str1,' kHz'])
subplot(2,1,2),plot(tfast,yfast,'g')
hold on
stem(t(1:11),x(1:11))
hold off
str2 = num2str(harmsamp);
title(['Resulting frequency ',str2,' kHz'])
sound(xhi,fhi)
sound(x,fs)

pause
end