% E04_7.m                26oct15   07apr21
% simple filters
clear all; close all
return;

%%
 num = [1 1 1 1 1 1 1 1 1]/9;    den = [1];           %  04.2-14
 num = [1  -1.2  1];    den = [1  -1.3  1.04  -.222]  %  04.2-16
% num = [1];    den = [1  -1.84  .85];   % stable 
% num = [1];    den = [1  -1.85  .85];   % unstable,  04.2-37
% num = [1];    den = [1  -2  1];        % unstable,  
 fvtool(num,den);

%%  what if the delay is not an integer?    plot for slide 04.3-47
n = linspace(0,10,100); 
s = 2*n .* exp(-.8*n); 
subplot(311); plot(n,s); hold on; stem(n(1:10:end),s(1:10:end)); hold off;
s1 = [zeros(1,10) s(1:end-10)];
subplot(312); plot(n,s1); hold on; stem(n(1:10:end),s1(1:10:end)); hold off;
s2 = [zeros(1,6) s(1:end-6)];
subplot(313); plot(n,s2); hold on; stem(n(1:10:end),s2(1:10:end)); hold off;

%%  truncated-sinc FIR impulse response
N = 34;    % e.g. range [20,200]   %  04.3-54 
n = [-N:N];
omc = .33 * pi;
h = sin(omc*n) ./ (pi*n);  
h(N+1) = omc/pi;   % avoid "NaN" in n=0
H = fft(h);  Hm = abs(H);  Hp = unwrap(angle(H));
subplot(311); stem(n,h); axis tight;
subplot(312); stem([0:N]/N*pi,Hm(1:N+1)); axis tight;
subplot(313); stem([0:N]/N*pi,Hp(1:N+1)); axis tight;
          % note linear phase in the passband
 fvtool(h,1);

%% moving average and modified symmetric filter
b1 = [1  1  1  1  1  1  1]/7;     a1 = [1];      %  04.5-22
b2 = [.5  1  1  1  1  1  .5]/6;    a2 = [1];
 fvtool(b2,a2,b1,a1);

%% simple FIR lowpass filter, and its cascaded version
num1 = [1 1]/2;  den1 = 1;        %  04.6-12
ze = [-1 -1 -1 -1]';   po = 0;    %  4 cascaded stages: 4 zeros in -1
[num,den] = zp2tf(ze,po,1);  
num = num ./ sum(num);            %  normalized filter coeffs.
 fvtool(num1,den1); fvtool(num,den);          %  compare

%% switch from lowpass to highpass filter
num1 = [1  4  6  4  1];   %  lowpass filter  04.6-17
num2 = [1  -4  6  -4  1];   %   changed to highpass
   fvtool(num1,1,num2,1);    %  compare
%%
 num = [1 0 -1];  fvtool(num,1);  % Central-difference differentiator (04.6-55)
	
%% simple IIR filters
%%%%%%% homework: y(n) = (1-a)*x(n) + a*y(n-1),  a = [-1,1)
%
 for alpha = -.9:.1:.9;   beta = .31;
% for beta = -.9:.1:.9;   alpha = .52;
  b = [1  1] * (1-alpha)/2;    a = [1  -alpha];    % lowpass, 04.6-27
% b = [1  -1] * (1+alpha)/2;    a = [1  -alpha];   % highpass, 04.6-231
% b = [1  0  -1] * (1-alpha)/2;   a = [1  -beta*(1+alpha)  alpha];   % bandpass, 04.6-39
% b = [1  -2*beta  1]* (1+alpha)/2;   a = [1  -beta*(1+alpha)   alpha];  % bandstop, 04.6-43
   subplot(222);  zplane(b,a);    % zeros and poles
   ome = [0:.01:pi];  H = freqz(b,a,ome); Hm = abs(H); Hp = angle(H);
   subplot(221);  plot(ome,Hm);  axis tight;    %  magnitude of freq. response
   subplot(223);  plot(ome,Hp);  axis tight;    %  phase of freq. response
   subplot(224);  stem(impz(b,a,30));  axis tight;    %  impulse response
 pause()
 end

%% example; effects of filter bandwidth
n = [0:299];
x =  sin(.1*pi*n) + cos(.17*pi*n) + randn(size(n))*10;
 b1 = [1 1]/2; a1 = 1;   b2 = [1 4 6 4 1]/16; a2 = 1;
%   alpha = .45;   beta = cos(.1*pi);  % bandpass filter, centered on cos(.1pi)
%   b1 = [1  0  -1] * (1-alpha) /2;    a1 = [1  -beta*(1+alpha)  alpha];   %  04.6-32
%   alpha = .95;   beta = cos(.1*pi);  % bandpass filter, centered on cos(.1pi)
%   b2 = [1  0  -1] * (1-alpha) /2;    a2 = [1  -beta*(1+alpha)  alpha];   
y1 = x;  y2 = x;
for i=1:5    % cascaded filters (e.g. 5)
    y1 = filter(b1,a1,y1);
    y2 = filter(b2,a2,y2);
    end;
subplot(311); plot(x);
subplot(312); plot(y1);
subplot(313); plot(y2);

%% compare IIR with FIR filter
clf
 b1 = [1 2];    a1 = [1 .4 -.12];       % IIR
 b1 = [1 2];    a1 = [1 -.7 .6];        % IIR
L = 13;                                 % L = 3 ... 11 
b2 = impz(b1,a1,L)';   a2 = [1, zeros(1,L-1)];      % similar FIR, using long division   (see E03_x.m)
   % note better make num and den same order (L-1 poles in z=0)
   % note impz output is a column vector, needs transpose
subplot(231);  zplane(b1,a1);    % IIR zeros and poles
subplot(232);  zplane(b2,a2);    % FIR zeros
ome = [0:.01:pi];  H1m = abs(freqz(b1,a1,ome));  H2m = abs(freqz(b2,a2,ome));  
 subplot(234);  plot(ome,H1m,ome,H2m);  axis tight;  % freq. responses (magn.)
 H1p = angle(freqz(b1,a1,ome));  H2p = angle(freqz(b2,a2,ome));  
 subplot(235);  plot(ome,H1p,ome,H2p);  axis tight;  % freq. responses (phase)

%% and with linear phase FIR filter
[ze,po,k] = tf2zp(b2,a2);
ze_linphase = [ze; ones(L-1,1)./abs(ze).*exp(-j*angle(ze))]; 
po_linphase = [po; po];
[b3,a3] = zp2tf(ze_linphase,po_linphase,k);
subplot(231);  zplane(b1,a1);    % IIR zeros and poles
subplot(232);  zplane(b2,a2);    % FIR zeros
subplot(233);  zplane(b3,a3);    % linear phase FIR zeros
H3m = abs(freqz(b3,a3,ome));  H3p = angle(freqz(b3,a3,ome));  
subplot(234);  plot(ome,H1m,ome,H2m,ome,sqrt(H3m*prod(abs(ze))),'k');  axis tight;    %  freq. responses
 % note: "sqrt" and "prod(abs(ze))" compensate gain of H3 for the presence of the 1/r zeros
 % each mirrored zero makes the gain of the pair = (1-r)^2/r instead of (1-r) 
 subplot(235);  plot(ome,H1p,ome,H2p);  axis tight;  % freq. responses (magn.)
 subplot(236);  plot(ome,H1p,ome,H3p,'k',ome,unwrap(H3p),'k--');  axis tight;  % freq. responses (phase)
 
