% Create an all-pass filter with phase distortion
c = [-1.4, 0.9];
Hd_ap = dfilt.allpass(c);
%phasez(Hd_ap);

% Create a lowpass filter
Fpass = 0.45;            % Passband Frequency
Fstop = 0.5;             % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor
b  = firgr('minorder', [0 Fpass Fstop 1], [1 1 0 0], [Dpass Dstop], ...
           {dens});
Hd_lp = dfilt.dffir(b);

% Cascade filter
Hd_casc = cascade(Hd_lp, Hd_ap);

% Signal and allpass filtered signal
n = 1:10000;
nh = 4;  T = 10;
y = zeros(size(n));
for k=1:nh, y(k:T:end) = 1; end
out_ap = filter(Hd_ap, y);
fh = figure;
sp = 2000;  np = 100;
subplot(2,1,1);  stem(y(sp:sp+np),'filled');
title('Signal');
subplot(2,1,2);  stem(out_ap(sp:sp+np),'filled');
title('APF signal');
phasez(Hd_ap);
pause;

% Square wave signal with filter applied
n = 1:10000;
nh = 4;
fh2 = figure;
for T = [8 9 10 11 12]
  
  % Create signal
  y = zeros(size(n));
  for k=1:nh, y(k:T:end) = 1; end
  
  % Filter signal
  out_ap = filter(Hd_ap, y);
  out_lp = filter(Hd_lp, y);
  out_casc = filter(Hd_casc, y);

  % Plots
  sp = 2000;  np = 100;
  subplot(3,1,1);  stem(y(sp:sp+np),'filled');  
  title('Signal');
  subplot(3,1,2);  stem(out_lp(sp:sp+np),'filled');  
  title('Linear phase LPF');
  subplot(3,1,3);  stem(out_casc(sp:sp+np),'filled');  
  title('Nonlinear phase LPF');

  pause;
end
