% filters used for speech synthesis
% http://uvafon.hum.uva.nl/david/ma_ssp/2010/Klatt-1980-JAS000971.pdf
% impulse invariant design

clear all
Fs = 8000 ;
T = 1/Fs ;
tmax = 4;
t=0:T:tmax ;

out1 = zeros(size(t));
x = chirp(t, 1, tmax, Fs/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% resonator -- bandpass
% resonant amplitude out is approx F/BW
F = 500 ; % Hz
BW = 30 ; % Hz
C = -exp(-2*pi*BW*T);
B = 2 * exp(-pi*BW*T) * cos(2*pi*F*T) ;
A = 1 - B - C ;

% Use A, B, C to compute filter response
 for i=3:length(t)
     out1(i) = A*x(i) + B*out1(i-1) + C*out1(i-2) ;
 end

% plot freq response
figure(1); clf
subplot(2,1,1)
[h,w] = freqz(A, [1 -B -C], 1000, Fs);
loglog(w,abs(h),'b', 'linewidth', 2)
hold on
[h,w] = freqz(fix(256*A)/256, [1 -fix(256*B)/256 -fix(256*C)/256], 1000, Fs);
loglog(w,abs(h),'r', 'linewidth', 2)
legend('exact', '8:8')
title(['F=', num2str(F), '   BW=', num2str(BW)])
set(gca, 'xlim', [10 4000])
subplot(2,1,2)
plot(t,out1);
xlabel('time (freq of chirp in KHz)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% lopass 
%6 db cutoff is at BW/2
F = 0 ; % Hz
BW = 1000 ; % Hz cutoff=1000
C = -exp(-2*pi*BW*T);
B = 2 * exp(-pi*BW*T) * cos(2*pi*F*T) ;
A = 1 - B - C ;

% Use A, B, C to compute filter response
 for i=3:length(t)
     out1(i) = A*x(i) + B*out1(i-1) + C*out1(i-2) ;
 end
 
figure(2); clf
subplot(2,1,1)
[h,w] = freqz(A, [1 -B -C], 1000, Fs);
loglog(w,abs(h),'b', 'linewidth', 2)
hold on
[h,w] = freqz(fix(256*A)/256, [1 -fix(256*B)/256 -fix(256*C)/256], 1000, Fs);
loglog(w,abs(h),'r', 'linewidth', 2)
legend('exact', '8:8')
set(gca, 'xlim', [10 4000])
title(['F=', num2str(F), '   BW=', num2str(BW), '  6db cutoff=', num2str(BW/2)])
subplot(2,1,2)
plot(t,out1);
xlabel('time (freq of chirp in KHz)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


