% downsamp_movie.m
%
% Animation showing downsampled signal in time and frequency domains as
% the input signal frequency is swept.
%
% 2015/03/04  Added more information printed to screen and fixed errors
%             on some figure titles.
% 2012/03/13  Re-written. shows time domain waveforms now among other things
%             psd() gives annoying warning message so I switched to pwelch()
% 2008/03/03  Written


clear;          % needed when changing waveforms

Length          = 2^10;   % 2^18
NumFrames       = 64;     % 128

Ymax            = +1.1;
Ymin            = -1.1;
MagicScaling    = 0.55;

fprintf(1, '===== downsamp_movie.m\n');
fprintf(1, '  Shows two time-domain views and two frequency-domain views of\n');
fprintf(1, '  a downsampled-by-2 cosine as it is swept in frequency from\n');
fprintf(1, '  zero (DC) to the Nyquist frequency pi.\n');
fprintf(1, '  Shows wnsamp_move.m\n\n');
fprintf(1, '  Set the size and location of the four figures on your screen\n');
fprintf(1, '    so that they are easily visible.\n\n');
fprintf(1, '  Note the downsampled waveform goes to zero occasionally only\n');
fprintf(1, '    because of the highly-unusual case here of a perfect tone \n');
fprintf(1, '    and where the sequence length is always a perfect multiple\n');
fprintf(1, '    of the frequency.\n\n');
fprintf(1, '  Increase value of variable "NumFrames" for smoother animations.\n');
fprintf(1, '    Decrease it to run faster.\n\n');
fprintf(1, '  Length of signal sequences: %i\n', Length);
fprintf(1, '  Number of frames:           %i\n\n', NumFrames);

%----- Main loop
freq_step       = Length/NumFrames;
for freq = 0:freq_step:Length
   in           = cos(freq*1*pi*(1:Length)/Length);
   in_down2_odd = in(1:2:end);          % keeps odd samples
   in_down2_even= in(2:2:end);          % keeps even samples

   % Plot time-domain waveforms (first few samples)
   PlotSampleRange = 25;
   figure(1); clf;
   subplot(3, 1, 1), stem(0:PlotSampleRange-1, in(1:PlotSampleRange));
     title('Input waveform--first few samples');
     axis([0 (PlotSampleRange-1) Ymin Ymax]);
   subplot(3, 1, 2), stem(0:PlotSampleRange-1, in_down2_odd(1:PlotSampleRange));
     title('Input waveform downsampled by 2--first few Odd samples');
     axis([0 (PlotSampleRange-1) Ymin Ymax]);
   subplot(3, 1, 3), 
     stem(0:PlotSampleRange-1, in_down2_even(1:PlotSampleRange));
     title('Input waveform downsampled by 2--first few Even samples');
     axis([0 (PlotSampleRange-1) Ymin Ymax]);

   % Plot time-domain waveforms (all)
   figure(2); clf;
   subplot(3, 1, 1), plot(0:(Length-1), in);
     title('Input waveform--all samples');
     axis([0 (Length-1) Ymin Ymax]);
   subplot(3, 1, 2), plot(0:(Length/2-1), in_down2_odd);
     title('Input waveform downsampled by 2--all odd samples');
     axis([0 (Length-1) Ymin Ymax]);
   subplot(3, 1, 3), plot(0:(Length/2-1), in_down2_even);
     title('Input waveform downsampled by 2--all even samples');
     axis([0 (Length-1) Ymin Ymax]);

   % Plot frequency-domain spectrum using psd()
   figure(3); clf;
   subplot(3, 1, 1), pwelch(in);
     title('Frequency-domain spectrum using psd()--input [0 - pi]');
     axis([0 1 -50 20]);
   subplot(3, 1, 2), pwelch(in_down2_odd);
     title('Frequency-domain spectrum using psd()--Downsampled odd samples');
     axis([0 1 -50 20]);
     axis([0 1 -50 20]);
     axis([0 1 -50 20]);
   subplot(3, 1, 3), pwelch(in_down2_even);
     title('Frequency-domain spectrum using psd()--Downsampled even samples');
     axis([0 1 -50 20]);

   % Plot frequency-domain spectrum using abs(fft())
   figure(4); clf;
   subplot(3, 1, 1), plot(abs(fft(in)));
     title('Frequency-domain spectrum using abs(fft())--input [0 - 2pi]');
     axis([0 Length 0 Length*MagicScaling]);
   subplot(3, 1, 2), plot(abs(fft(in_down2_odd)));
     title('Frequency-domain spectrum using abs(fft())--Downsampled odd samples');
     axis([0 Length 0 Length*MagicScaling]);
   subplot(3, 1, 3), plot(abs(fft(in_down2_even)));
     title('Frequency-domain spectrum using abs(fft())--Downsampled even samples');
     axis([0 Length 0 Length*MagicScaling]);

   if (freq==0 ...
     || freq==freq_step ...
     || freq==2*freq_step ...
     || freq==3*freq_step ...
     || freq==Length-freq_step)
    %|| freq==(Length/2)-freq_step ...
    %|| freq==(Length/2) ...
    %|| freq==(Length/2)+freq_step ...
     fprintf(1,'Step = %i/%i; Freq = %i/%i,  ', freq/freq_step, ...
        (Length/freq_step -1), freq, Length);
     junk = input('--Press return to continue--');
     end
end

