% Demonstration of frequency response of system for unit amplitude cosine
% inputs at different frequencies

% Frequencies and time instants for plots
wv = -5:0.01:5;
tv = -5:0.1:5;
wcv = 0:0.5:5;

ncase = 2;  %% ----- change this for different filters -----
switch ncase

  case 0
    % First order lowpass filter
    R = 1;  C = 1;
    Hv = 1./(1+j*wv*R*C);
    
  case 1
    % Moving average filter
    tau = 6;
    Hv = sin(tau*wv/2)./(tau*(wv)/2);

  case 2
    % Ideal lowpass filter
    wcutoff = 2;
    td = 0.1;  % time delay (linear phase)
    Hv = exp(-j*wv.*td);  Hv(find(abs(wv)>wcutoff)) = 0;
  
end

% Plot frequency response magnitude and phase
fh = figure;
subplot(3,1,1);
plot(wv,abs(Hv));  axis([-5 5 0 1]);
xlabel('\omega');  ylabel('|H|');
subplot(3,1,2);
plot(wv,angle(Hv));  axis([-5 5 -2 2]);
xlabel('\omega');  ylabel('arg(H)');

% Demonstration for each cosine frequency
lhv = [];
for i=1:length(wcv)
  wc = wcv(i);  % current cosine frequency
  xv = cos(wc*tv);  % current cosine input signal
  Hwc = Hv(find(wv==wc));  % current system frequency response value (complex)
  yv = abs(Hwc)*cos(wc*tv + angle(Hwc));  % current output signal
  
  % Delete old plot handles
  if ~isempty(lhv)
    for hi=1:length(lhv), delete(lhv(hi)); end
    lhv = [];
  end
  
  % Updated plots
  subplot(3,1,1);  
  lh1 = line([wc; wc], [0; 1]);  set(lh1,'Color',[0 1 0]);
  lh2 = line([-wc; -wc], [0; 1]);  set(lh2,'Color',[0 1 0]);
  subplot(3,1,2);  
  lh3 = line([wc; wc], [-2; 2]);  set(lh3,'Color',[0 1 0]);
  lh4 = line([-wc; -wc], [-2; 2]);  set(lh4,'Color',[0 1 0]);
  lhb = line([-5; 5], [0; 0]);  set(lhb,'Color','k');
  lhv = [lh1 lh2 lh3 lh4];

  subplot(3,1,3);
  plot(tv,xv,'r',tv,yv,'g');  axis([-5 5 -1 1]);
  lha = line([-5; 5], [0; 0]);  set(lha,'Color','k');
  xlabel('t');  
  legend('Input x(t)','Output y(t)');
  pause;
end