% Matt Montag
% EEN 540 Project 2 
% Problem 1: Vocal Tract Area Function
close all;

A = [5 5 5 5 6.5 8 8 8 8 8 8 8 8 6.5 5 4 3.2 1.6 2.6 2.6 2 1.6 1.3 1 0.65 0.65 0.65 1 1.6 2.6 4 1 00.2 1.3 1.6 2.6];
E = [8 8 5 5 4 2.6 2 2.6 2.6 3.2 4 4 4 5 5 6.5 8 6.5 8 10.5 10.5 10.5 10.5 10.5 8 8 6.5 6.5 6.5 6.5 1.3 1.6 2 2.6];
I = [4 4 3.2 1.6 1.3 1 0.65 0.65 0.65 0.65 0.65 0.65 0.65 1.3 2.6 4 6.5 8 8 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 8 8 2 2 2.6 3.2];
O = [3.2 3.2 3.2 3.2 6 6 6 6 6 6 6 6.5 13 13 16 13 10.5 10.5 8 8 6.5 6.5 5 5 4 3.2 2 1.6 2.6 1.3 0.65 0.65 1 1 1.3 1.6 2 3.2 4 5 5 1.3 1.3 1.6 2.6];
U = [0.65 0.65 00.4 00.4 00.4 0.6 2 5 7 9 11 13 13 10.5 8 6.5 5 3.2 2.6 2 2 2 1.6 1.3 2 1.6 1 1 1 1.3 1.6 3.2 5 8 8 10.5 10.5 10.5 2 2 2.6 2.6];
vowels = { A E I O U ; 'A' 'E' 'I' 'O' 'U' };
voices = { 'male' 'female' 'child'; 120 240 300; 1 0.8 0.6; .85 .7 .6; .2 .3 .9 }; % label; pitch; vocal tract scale; glottal tense; breathiness
tube_width = 0.5; % (cm)
fs = 44100;
dur = 1; % duration of synthesized sample (sec)
c = 345; % speed of sound (m/s)
denlossy = {};
% Tube modeling
for voice=1:3 % number of different voices
	person = voices{1,voice}; % a label for plots and filenames
	pitch = voices{2,voice}; % (hz)
	vocal_tract_scale = voices{3,voice}; % as a fraction of male vocal tract length
	max_length = 17.5 * vocal_tract_scale; 
	glottal_tense = voices{4,voice};
	% Rosenberg pulse
	signal = [];
	for i=1:6
		signal = [signal rosenberg(glottal_tense, .65, pitch, fs)];
	end
	figure; subplot(2,1,1); plot(signal); ylim([-1 2]); xlabel('Samples','fontsize',7);
	title('General glottal excitation using Rosenberg pulse');
	Y = fft(signal);
	Y = abs(Y(1:floor(length(Y)/2)));
	F = linspace(0,fs/2,length(Y))';
	subplot(2,1,2);
	semilogy(F,Y); ylabel('dB'); xlabel('Frequency','fontsize',7); xlim([0 fs/2]);
	saveas(gcf,[person '-glottalpulse'],'png');
	for i=1:length(vowels)
		vow = fliplr(vowels{1,i});
		% Interpolate tubes based on vocal tract length and bandwidth
		N = round(vocal_tract_scale * length(vow) * tube_width / 100 *  fs/2 * 4 / c);
		vow = interp1(vow, linspace(1, length(vow), N));
		vow = [vow 25 25]; % append "dummy tubes" for display only
		subplot(2,2,1);
		stairs(vow); ylim([-26 26]); hold on;
		stairs(-vow); hold off;
		title(['"' vowels{2,i} '" - Vocal tract area function']);
		% Calculate reflection coefficients
		refs = []; % calculate reflection coefficients for each boundary
		for j=1:length(vow)-3 %omit the last 2 "dummy tubes" 
			% r = (next-this) / (next+this)
			ref = (vow(j+1)-vow(j)) / (vow(j+1)+vow(j));
			refs = [refs ref];
		end
		refs = [refs 1]; % add the infinite tube
		refslossy = [refs(1:end-1) .71]; % add a finite (lossy) tube
		subplot(2,2,2);
		stem(refs); ylim([-1 1]);
		title(['"' vowels{2,i} '" - Reflection coefficients']);
		% Durbin's recursion: from reflection coeffs to filter coefficients
		den = [1]; denlossy{i} = [1];
		for k=1:length(refs)
			den = [den 0] + [0 fliplr(refs(k)*den)];
			denlossy{i} = [denlossy{i} 0] + [0 fliplr(refslossy(k)*denlossy{i})];
		end
		num = [zeros(1,round(length(refs)/2)) 1];
		subplot(2,1,2);
		[h,f] = freqz(num,den,512,fs);
		plot(f,20*log10(abs(h))); xlim([0 fs/2]);
		hold on; 
		[h,f] = freqz(num,denlossy{i},512,fs );
		plot(f,20*log10(abs(h)),'r'); xlabel('Frequency','fontsize',7);
		title(['"' vowels{2,i} '" - Magnitude response']);
		% Radiation model
		alpha = .8; R = [1 -alpha]; num = conv(R, num);
		[h,f] = freqz(num,denlossy{i},512,fs );
		plot(f,20*log10(abs(h)),'m'); hold off; legend('Lossless','Lossy','Lossy+Radiation');
		saveas(gcf,[person '-' vowels{2,i}],'png');
		%break;
		pause(0.1);
	end
	% Synthesized vowel sounds
	signal = [];
	
	for i=5.5:-1:1
		signal = [signal .6*rosenberg(.9, .4, pitch, fs) zeros(1,2^i+rand*fs/200)];
	end
% 	for i=1:pitch * dur
% 		signal = [signal rosenberg(glottal_tense, .6, pitch, fs)];
% 	end
	pitchdip = (.98+rand*.02);
	pvec = [linspace(pitch,pitch*pitchdip,pitch*dur*14/16) linspace(pitch*pitchdip,pitch/4,pitch*dur*1/10)];
	avec = pvec.*pvec/pitch;
	gvec = glottal_tense*pvec/pitch;
	for i=1:length(pvec)
		roz = avec(i)*rosenberg(gvec(i), .6, pvec(i), fs);
		noisy = conv(rand(1,length(roz)-15), hamming(16));
		roz = (1+noisy*0.03*voices{5,voice}) .* roz;
		signal = [signal roz];
	end
	env = [linspace(0,1,length(signal)/16) ones(1,length(signal)*.75) linspace(1,0,length(signal)/8)];
	signal = conv(signal, [.6 1 1 .6]);
	[env, signal] = pad(env, signal);
	signal = (signal-.5) .* env;
	for i=1:length(denlossy)
		den = denlossy{i};
		output = filter(num, den, signal);
		output = output .* 0.7 / max(abs(output)); %normalize
		wavwrite(output, fs, 16, [person '-' vowels{2,i} '.wav']);
		soundsc(output, fs);
	end
end

%%
% 
% % test a sliding filter.
% signal=[];
% for i=1:100
% 		signal = [signal rosenberg(glottal_tense, .65, pitch, fs)];
% end
% steps = 1000;
% output = zeros(1,40000);
% for i=1:steps
% 	AA = i/100;
% 	BB = 1-AA;
% 	den = BB*denlossy{3} + AA*denlossy{4};
% end