%  This script graphically illustrates aliasing.  It creates/plays a movie-like plot
%  showing a densly sampled (non-aliased) sinusoid with its subsampled aliased version
%  overlayed with a stem plot.  The sampling is applied over a range of frequencies
%  can can be selected along with the sampling rates in the first several lines of the
%  script.  The script stores frames of the figure in a movie file and saves it as an avi
%  file in the current directory for replay with any avi player
%
% Written by Kevin D. Donohue (donohue@engr.uky.edu) updated January 2009

fs1 = 200;  % Sampling Frequency to produce aliasing
fs2 = 4000; % Sampling with no aliasing (must be at least greater than 2*max(f1)
            %  best if 10 to 100 times greater)
fbegin = 350;  % Initial Signal Frequency
fend = 450;  % Final Signal Frequency
finc = 2;  % Frequency sequence increment
int = .1;   % Time of interval to sample the plotted signal in seconds (should
            % be less than one for good graphic quality).
t1 = [0:fix(fs1*int)]/fs1;  %  Graphing axis for fs1 Hz sampling
t2 = [0:fix(fs2*int)]/fs2;  %  Graphing axis for fs2 Hz sampling

mindex = 0;  % Initalize index for movie file
%  Loop indexes are frequencies ploted to show alising
h1 = figure;  %  Set up figure for frames to display
set(h1,'Position', [137 299 627 385])
set(h1,'NextPlot', 'Replace' ) ;   %  Next time replace with new axis info
%  Step through frequency range, sample each, and create movie frame
for k=fbegin:finc:fend   
	f1 = k;   % Frequency of sampled waveform
	s1 = cos(2*pi*f1*t1);  %  Aliased signal
	s2 = cos(2*pi*f1*t2);  %  Signal with no aliasing (if 2*f1 < fs2) 
	%  Compute aliased frequency
    fal = mod(f1,fs1);  %  Remove multiple of sampling rate
	if fal >= fs1/2   %  if greater than half the sampling frequency fold back
        fal = fs1-fal;
	end
    s3 = cos(2*pi*fal*t2);  %  Over sampled aliased signal for plotting
    %  Plot to illustrate the aliasing 
	p1 = stem(t1,s1,'b');  %  Plot aliased signal with stem plot
    ah1 = get(h1,'CurrentAxes');
	set(ah1, 'FontSize', 12, 'YLim', [-1.1 1.1],'NextPlot', 'Add' ) ;
    plot(t2,s3,'b-','Linewidth',3);  %Overlay outline of aliased signal  
	p2 = plot(t2,s2,'r--');  %  Overlay non-aliased signal
	set(ah1,'NextPlot', 'Replace' ) ;   %  Next time replace with new axis info
    set(p1, 'LineWidth', 3)
	set(p2, 'LineWidth', 1)
	xlabel('Seconds','FontSize', 14)
	ylabel('Amplitude','FontSize', 14)
	    
	title(['Red: ' num2str(f1) ' Hz, Blue Aliased to: ' num2str(fal) ' Hz '],'FontSize', 14)
	
    mindex=mindex+1;   %  increment index for movie file
    mv(mindex) = getframe(h1);  % Movie file, use movie2avi(mv) to save as avi file
    pause(.2)   %  Pause in loop so you can see the graphic
end
movie2avi(mv,'alias2.avi','fps',2,'compression','none');