% Demonstration of aliasing (SAGE 2000 -- L. Braile)
% generate 2 sinusoids
dt=0.01;  % 100 samples per second, Nyquist = 50 HZ
t=[0:dt:1.0];
f1=12;  % 12 HZ sinusoid
s1=sin(2*pi*t*f1);
f2=30;  % 30 HZ sinusoid
s2=sin(2*pi*t*f2);
s3=s1+s2;
%plot sinusoids
hold on
plot(t,s1,'-r','linewidth',2)
plot(t,s2,'--b','linewidth',2)
plot(t,s3,':k','linewidth',2)
hold off
set(gca,'fontsize',16,'linewidth',2)
xlabel('Time (s)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('12 Hz and 30 Hz sinusoids and sum, sampled at 0.01 s')
figure
% calculate and plot 1024 point fourier transform
S = fft(s3,1024);
Sxx = S.*conj(S)/1024; % calculate periodogram
f=100*(0:512)/1024;  % calculate frequencies (513) contained in 
                     % independent part of 1024 point spectrum Sxx.
plot(f,Sxx(1:513))   % plot spectrum (513 points from 0 to fnyq Hz)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
set(gca,'fontsize',16,'linewidth',2)
title('FFT of sum of 12 and 30 Hz sinusoids sampled at 0.01 s')
figure
% plot with log amplitude scale
semilogy(f,Sxx(1:513))
axis([0 50 .0001 10]);
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
set(gca,'fontsize',16,'linewidth',2)
title('FFT of sum of 12 and 30 Hz sinusoids sampled at 0.01 s')
figure
% Resample s3 with an incorrect sample interval
dt2 = 0.02; % 50 samples per second, Nyquist = 25
s4=s3(1:2:101);
t4=[0:0.02:1.0];
plot(t4,s4,'-r','linewidth',2)
xlabel('Time (s)','fontsize',16)
ylabel('Amplitude','fontsize',16)
set(gca,'fontsize',16,'linewidth',2)
title('Sum of 12 and 30 Hz sinusoids re-sampled at 0.02 s')
figure
% calculate and plot 1024 point fourier transform
S = fft(s4,1024);
Sxx4 = S.*conj(S)/1024;
f=50*(0:512)/1024;
plot(f,Sxx4(1:513))
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
set(gca,'fontsize',16,'linewidth',2)
title('FFT of sum of 12 and 30 Hz sinusoids resampled at 0.02 s')
figure
% plot with log amplitude scale
semilogy(f,Sxx4(1:513))
axis([0 25 .0001 10]);
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
set(gca,'fontsize',16,'linewidth',2)
title('FFT of sum of 12 and 30 Hz sinusoids resampled at 0.02 s')