%speech

% 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;

% since BW*T and F*T appear only as products
% normalize to reduce dimensionality
F = 0.02:0.02:1/2 ; % F/(Fs/2)
BW = F/10 ; % set bandwidth as a fraction of F
T = 1;

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

% approximate Cm B and A by polynomials
pc = polyfit(F, C, 2); % works within 2e-4
pb = polyfit(F, B, 3); % works within 0.02; within 0.01 for F<0.4
Cp = polyval(pc, F); ;
Bp = polyval(pb, F);
Ap = 1-polyval(pc, F)-polyval(pb, F);

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

% plot errors between A, B and C and approximations
figure(2);clf;
hold on
plot(F, C-Cp, 'bx')
plot(F, B-Bp, 'rx')
plot(F, 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(F)
    [h,w] = freqz(A(i), [1 -B(i) -C(i)], 500,1);
    loglog(w,abs(h),'b', 'linewidth', 1)
    hold on
    [h,w] = freqz(Ap(i), [1 -Bp(i) -Cp(i)], 500,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], 500,1);
    loglog(w,abs(h),'g', 'linewidth', 1)

end