%LOW PASS second order impluse invariant design

% search: "speech formant synthesis  algorithm"
% http://uvafon.hum.uva.nl/david/ma_ssp/2010/Klatt-1980-JAS000971.pdf
% http://en.wikipedia.org/wiki/Speech_synthesis

% Fs = 8000;
% T = 1/Fs;
% F = 0:0.05:Fs/2 ;
% BW = F/50;

clear all
% since BW*T and F*T appear only as products
% normalize to reduce dimensionality
BW = 0:0.05:1/2; % set bandwidth for lowpass
T = 1;

C = -exp(-2*pi*BW*T);
B = 2 * exp(-pi*BW*T)  ;
A = 1 - B - C ;

% approximate Cm B and A by polynomials
pc = polyfit(BW, C, 3); % 
pb = polyfit(BW, B, 3); % 
Cp = polyval(pc, BW); ;
Bp = polyval(pb, BW);
Ap = 1-polyval(pc, BW)-polyval(pb, BW);

% compare A, B and C with approximations
figure(1);clf;
plot(BW, C, 'b')
hold on
plot(BW, Cp, 'bx')
plot(BW, B, 'r')
plot(BW, Bp, 'rx')
plot(BW, A, 'g')
plot(BW, Ap, 'gx')
legend ('C', 'B', 'A')

% plot errors between A, B and C and approximations
figure(2);clf;
hold on
plot(BW, C-Cp, 'bx')
plot(BW, B-Bp, 'rx')
plot(BW, A-Ap, 'gx')
legend ('C', 'B', 'A')

% plot frequency response for exact, approximated, and 8:8 fixed truncation
figure(3); clf; 
for i=1:length(BW)
    [h,w] = freqz(A(i), [1 -B(i) -C(i)], 1000,1);
    loglog(w,abs(h),'b', 'linewidth', 1)
    hold on
    [h,w] = freqz(Ap(i), [1 -Bp(i) -Cp(i)], 1000,1);
    loglog(w,abs(h),'r', 'linewidth', 1)
    [h,w] = freqz(fix(256*Ap(i))/256, [1 -fix(256*Bp(i))/256 -fix(256*Cp(i))/256], 1000,1);
    loglog(w,abs(h),'g', 'linewidth', 1)

end