% Script that dynamically illustrates aliasing by creating an animation of a
% sampled sinewave with variable samplng rate, starting well above the
% Nyquist frequency and ending well below.

SignalFrequency=100; % Frequency of the sinewave
clf

% Quickly vary the sampling rate from 600 to 300
for sf=600:-5:300 % SamplingFrequency
    dt=1/sf; % time interval between samples
    t = 0:dt:.2; % time values at which the waveform is sampled
    signal = sin(2*pi*SignalFrequency*t); % Sampled signal
    plot(t,signal,'.-');
    axis([0 max(t) -1 1])
    xlabel('time')
    ylabel('Signal amplitude')
    PeakCount=round(length(allpeaks(t,signal)));
    title(['Signal frequency = ' num2str(SignalFrequency) '     Sampling frequency = ' num2str(sf)   '    PeakCount=' num2str(PeakCount)])
    drawnow
    pause(.1)
end

% Slowly vary the sampling rate from 300 to 130
for sf=300:-1:130 % SamplingFrequency
    dt=1/sf; % time interval between samples
    t = dt:dt:.2; % time values at which the waveform is sampled
    signal = sin(2*pi*SignalFrequency*t); % Sampled signal
    plot(t,signal,'.-');
    axis([0 max(t) -1 1])
    xlabel('time')
    ylabel('Signal amplitude')
    PeakCount=round(length(allpeaks(t,signal)));
    title(['Signal frequency = ' num2str(SignalFrequency) '     Sampling frequency = ' num2str(sf)   '     Number of peaks = ' num2str(PeakCount)])
    drawnow
    pause(.1)
end

function P=allpeaks(x,y)
% Looks for data points in x,y that have a y value greater than the
% adjacent points, then prints out the peak number, x value, and y value of
% those points. (Noisy data must be smoothed first).
% Example:
%  x=2:38
%  y=abs(randn(size(x)));
%  P=allpeaks(x,y);
%  disp('           peak      position      height')
%  disp(P)
%  plot(x,y)
%  text(P(:,2),P(:,3),num2str(P(:,1))) % Numbers peaks
%
peak=0;
for k = 2:length(x)-1
    if y(k) > y(k-1)
        if y(k) > y(k+1)
            peak = peak + 1;
            P(peak,:)=[peak x(k) y(k)];
        end
    end
end
end