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

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

%  -- bandpass
% resonant amplitude out is approx F/BW
F = logspace(2.0,3.3,25)/(Fs/2) ; % normalized freq
BW = 0.15 ; % of freq
for i=1:length(F)
    [b{i}, a{i}] = butter(1,[F(i)-F(i)*(BW/2), F(i)+F(i)*(BW/2)] );
end
two14 = 2^14;

% plot freq response
figure(1); clf
for i=1:length(F)
    %subplot(2,1,1)
    [h,w] = freqz(b{i}, a{i}, 3000, Fs);
    loglog(w,abs(h),'b', 'linewidth',1)
    hold on
    [h,w] = freqz(fix(two14*b{i})/two14, fix(two14*a{i})/two14, 5000, Fs);
    loglog(w,abs(h),'r', 'linewidth', 1)
    legend('exact', '2:14', 'Location', 'East')
    %title(['F=', num2str(F), '   BW=', num2str(BW)])
    set(gca, 'xlim', [10 4000])
    set(gca, 'ylim', [.01 1])
    xlabel('frequency (Hz)')
    ylabel('log amplitude')
%      subplot(2,1,2)
%      plot(t,out1);
%      xlabel('time (freq of chirp in KHz)')
end
return
%make a textfile with GCC source code in it.
% for A, B, C arrays
float_2_fix = 2^14 ;
fname='Filter_bandpass_14bit.h';
fid = fopen(fname,'w');
fprintf(fid,'#define Filter_bank_size %d\r',length(F));
fprintf(fid,'const int b1[%d]={\r',length(F));
for i=1:length(F)-1
    fprintf(fid,' %7d,\r',fix(b{i}(1)*float_2_fix));
end
fprintf(fid,' %7d};\r',fix(b{end}(1)*float_2_fix));

fprintf(fid,'const int a2[%d]={\r',length(F));
for i=1:length(F)-1
    fprintf(fid,' %7d,\r',fix(a{i}(2)*float_2_fix));
end
fprintf(fid,' %7d};\r',fix(a{end}(2)*float_2_fix));

fprintf(fid,'const int a3[%d]={\r',length(F));
for i=1:length(F)-1
    fprintf(fid,' %7d,\r',fix(a{i}(3)*float_2_fix));
end
fprintf(fid,' %7d};\r',fix(a{end}(3)*float_2_fix));
fclose(fid);

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


