%
%  this program shows the effect of sampling a signal (using impulse sampling)
%  and then filtering the signal to recover the original signal.
%
%  The sampling is modelled by modulation.
%
%  x(t)-->[sample/demodulate]-->xs(t) --> [LP Filter] --> xr(t)
%
function dummy
%
%  determine the time/frequency scale
%
   kHz = 1000.0;
%   
%  determine the simulation sampling frequency (in Hz)
%
   Fs = 5000*kHz;  % kHz
%
%  set the time interval for the original function
%
   T = 2.0/kHz;    % seconds
   t = [-T/2:1/Fs:T/2];
%
%  determine the function to sample
%
  x = @(t) 1*(abs(t)<2/kHz)+0*(abs(t)>=2/kHz);    % rect function
  x = @(t) sinc(t/(0.2/kHz)).*sinc(t/(0.2/kHz));  % sinc^2 function
%
%  determine the effecitive sampling rate and frequency
%
  fo = 20.0*kHz;   % sampling frequency
%
%  determine the number of times to replicate the signal
%
  Nr = 5;
%
%  determine the maximum frequency to plot
%
  fmax = (Nr-0.5)*fo;
%
%  modulate the signal (fake sampling)
%
  xs = x(t);
%  
  for k = 1:Nr
     xs = xs + 2*x(t).*cos(2*pi*k*fo*t);
  end;
%
%  examine the signal in the frequency domain
%
  [xSmag,xSphase,xsdf,xsfrq] = ctft_fft(xs,1/Fs);
%
%  scale the way impulse sampling scales
%
  xSmag = xSmag*fo;
%
%  determine the lowpass filter
%
  Flp = 2*kHz;  % cutoff frequency in Hz
  [b,a] = butter(5,Flp/(Fs/2));
%
%  look at the frequency response of the lowpass filter
%
   F = linspace(0.1,fmax,1000);
   H = (1/fo)*freqz(b,a,F,Fs);
   H_dB = 20*log10(abs(H));
   plot_H_dB = [fliplr(H_dB) H_dB];
   plot_F = [-fliplr(F) F];
%  
%  filter the output (run the output through the filter)
%  
  xr = filtfilt(b,a,xs);
%
  figure;
  orient tall
  
  subplot(4,1,1);  plot(t*kHz,x(t),'-b');  grid; 
  xlabel('Time (ms)');
  ylabel('x(t)');
  title(['Sampling Rate f_s = ', num2str(fo/kHz), ' kHz, Lowpass Cutoff Frequency f_c = ', num2str(Flp/kHz),' kHz']);
    
  subplot(4,1,2); plot(xsfrq/kHz,xSmag); axis([-fmax/kHz fmax/kHz 0 max(xSmag)]); grid; 
  xlabel('Frequency (kHz)');
  ylabel('|X_s(f)|');

  subplot(4,1,3); plot( plot_F/kHz, plot_H_dB ); 
  axis([-fmax/kHz fmax/kHz min(plot_H_dB) max(plot_H_dB)]); grid; 
  xlabel('Frequency (kHz)'); 
  ylabel('|H(f)| dB');
 
  subplot(4,1,4); plot(t*kHz, xr, '-r', t*kHz, x(t), '--b'); grid;
  legend('x_r(t)','x(t)');
  xlabel('Time (ms)');
  ylabel('Final');
  
   figure; plot(t*kHz, xr, '-r', t*kHz, x(t), '--b'); grid;
  legend('x_r(t)','x(t)');
  xlabel('Time (ms)');
  ylabel('Final');
  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  function [Xmag,Xphase,df,frq] = ctft_fft(x,dt)
%
% Computes FFT of input vector, returning the magnitude and phase spectra and
% a frequency vector for display.  Includes arranging the frequencies properly, proper
% scaling and units, and unwrapping the phase vector.
%  BAF  revised fall 2004
%
% usage:  [XMag,XPhase,df,frq]=ctft_fft(x,dt)
%       x = function for which FFT is calculated
%       dt = time spacing bewteen samples of the function
%       XMag = vector containing magnitude of FFT (two-sided magnitude spectrum)
%       XPhase = vector containing phase of FFT (two-sided phase spectrum)
%       df = frequency spacing between FFT elements
%       frq = vector containing frequency values - for plotting
%
N=length(x);
df=1/N/dt;
% compute FFT, with proper scale factor
X1=dt*fft(x);
% rearrange for proper sequencing of neg/pos frequency components
X=fftshift(X1);
% compute mag and phase vector, with phase unwrapping
Xmag=abs(X);
Xphase=unwrap(angle(X));
% generate frequency vector associated with FFT vectors
frq=linspace((-N/2)*df, (N/2-1)*df, N);
  