% eda07_06
% use z-transforms to:
%    check if a filter is minumum phase
%    represent it as a cascade of length-2 filters
%    compute the inverse filter

% i used as imaginary unit
clear i;

% define a filter
g = [6, 5.5, 5.3, 5, 4, 3, 2]';
N = length(g);
gN = g(N);

% find roots of g
r = roots(flipud(g));

% check size of roots
badroots=0;
for j = [1:N-1]
    if( abs(r(j)) <= 1 )
        badroots=badroots+1;
    end
end
disp(sprintf('%d roots not outside the unit circle', badroots));

% rebuild filter from its roots
gp = zeros(1,1);
gp(1)=gN;
for j = [1:N-1]
    gp = conv(gp,[-r(j),1]');
end
% there should be no imiginary part
gp = real(gp);

% check that it worked
E1 = (g-gp)'*(g-gp);
disp(sprintf('error of filter rebuilt from its roots: %f',E1));

% construct inverse filter, gi, of length Ni
Ni = 50;
gi = zeros(Ni,1);
gi(1)=1/gN;
% filter cascade, one filter per loop
for j = [1:N-1]
    % construct inverse filter of a length-2 filter
    tmp = zeros(Ni,1);
    tmp(1) = 1;
    for k = [2:Ni]
        tmp(k)=tmp(k-1)/r(j);
    end
    tmp = -tmp/r(j);
    gi = conv(gi,tmp);
    gi=gi(1:Ni);
end
% delete imaginary part (which should be zero)
gi = real(gi);

% test quality of inverse filter
gig=conv(g,gi);
gig=gig(1:Ni);
d = zeros(Ni,1);
d(1)=1;
E2 = (gig-d)'*(gig-d);
disp(sprintf('error of inverse filter: %f',E2));

% plot results
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, min(g), max(g)] );
plot( [1:Ni]', [g',zeros(Ni-N,1)']', 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('g(j)');
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, min(gi), max(gi)] );
plot( [1:Ni]', gi, 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('ginv(j)');
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, -1.2, 1.2] );
plot( [1:Ni]', gig, 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('[g*ginv](j)');







    
