% Matt Montag
% EEN 540 Project 3 
% Linear Predictive Coding

close all; clear all; clc;
m2f = 1;
stylized = 0;
noiseGain = 1;

if(m2f == 0) % female to male
	filename = 'female_speech.wav';
	[wav, fs] = wavread(filename);
	N = round(20 * fs/160); %floor(0.02 * fs); % window length (samples)
	child = 0;
	if(stylized ~= 0) noiseGain = 4 * fs/16000; end
	order = (fs/1500) + 2 + 4;
	pitchNum = 1;
	pitchDen = 1;
	poleFactor = 1;
	maxpitch = 250;
	tukeyshape = 1; % 1 = hanning; 0 = rect
else 
	filename = 'male_speech.wav';
	[wav, fs] = wavread(filename);
	N =	round(20 * fs/160);
	child = 1;
	if(stylized ~= 0) noiseGain = 2; end
	order = (fs/1000) + 2 + 4;
	pitchNum = 1;
	pitchDen = 1;
	poleFactor = 1.5;
	maxpitch = 280;
	tukeyshape = 1; % 1 = hanning; 0 = rect
end
voicedThresh = 0.85; % Ratio of 1st and 2nd xcorr values to consider voiced
speed = pitchNum/pitchDen;
hop = floor(N *.5); %(1-tukeyshape/2));
N = floor(N*speed);
I = 0;

%Spectral separation of sines+noise. Add the noise back in at the end.
%[Sines, Noise] = sinesNoise(wav, fs);

output = zeros(round(length(wav) * 1.2),1);
errOutput = output;
iterator = 0;
haltOnIterator = 28;
halt = 0;

for i=1:hop:length(wav)-N
	iterator = iterator + 1;
	if i+N-1 > length(wav) 
		continue;
	end
	xx = tukeywin(N,tukeyshape) .* wav(i:i+N-1);

	axx=xcorr(xx,xx);  % Autocorrelation of frame
	[ma, W] = max(axx); % Max val/location
	mu(iterator) = axx(W+1)/axx(W); % Pre-Emphasis Factor
	tonality = mu(iterator);
	
	if(tonality > voicedThresh) voiced = 1; else voiced = 0; end;
	
	[a, pe] = lpc(xx, order);
	err = filter(a, 1, xx);

	maxpeak = max(abs(err));
% 	[peakX, peakY] = peakseek(errNew./hamming(length(errNew)),round(fs/maxpitch),0);	
	[peakY, peakX] = findpeaks(abs(wav(i:i+N-1)), ...
		'minpeakdistance',round(fs/maxpitch/(1+stylized)),'minpeakheight',maxpeak*.2); % Important tweaky factor, minpeakheight

	if(pitchDen/pitchNum ~= 1)
		%resample to stretch the signal out
		%errNew = resample(err, pitchDen, pitchNum);
		padd = N;
		errNew = zeros(padd+round(length(err)*pitchDen/pitchNum),1);
		err = [zeros(padd/2,1); err; zeros(padd/2,1)]; % pad original signal
		%then replace peak areas with peak areas from original signal
		for pk=1:length(peakX)-1 
			period = peakX(pk+1) - peakX(pk);
			reach = round(period);
			origX = peakX(pk);
			scaledX = round(origX * pitchDen / pitchNum);
			errNew(padd/2+scaledX-reach:padd/2+scaledX+reach) = hamming(reach*2+1).* err(padd/2+origX-reach:padd/2+origX+reach);
		end
		%chop off the padding and apply a small window to the entire synthetic signal
		errNew = errNew(padd/2+1:length(errNew)-padd/2);
		err = err(padd/2+1:length(err)-padd/2);
		errNew = tukeywin(length(errNew), .05).*errNew;
	else
		errNew = err;
	end
	
	%insert complementary impulses in error signal to increase pitch.
	if(m2f==1) 
		if(child == 0)
			for pk=1:length(peakX)-1
				period = peakX(pk+1) - peakX(pk);
				quarter = round(period/4);
				midpt = round((peakX(pk)+peakX(pk+1))/2);
				if (midpt - quarter) < 1 || peakX(pk)-quarter < 1 || peakX(pk+1)+quarter > length(errNew)
					continue;
				end
				errNew(midpt-quarter:midpt+quarter) = ...
					( errNew(peakX(pk)-quarter:peakX(pk)+quarter) + ...
					errNew(peakX(pk+1)-quarter:peakX(pk+1)+quarter) ) ...
					/2;	
			end
		else
			for pk=1:length(peakX)-1
				period = peakX(pk+1) - peakX(pk);
				quarter = round(period/6);
				pt13 = round((2*peakX(pk)+peakX(pk+1))/3);
				pt23 = round((peakX(pk)+2*peakX(pk+1))/3);
				if (peakX(pk)-quarter < 1 || peakX(pk+1)+quarter > length(errNew))
					continue;
				end
				errNew(pt13-quarter:pt13+quarter) = ...
					( 2*errNew(peakX(pk)-quarter:peakX(pk)+quarter) + ...
					errNew(peakX(pk+1)-quarter:peakX(pk+1)+quarter) ) ...
					/3;		
				errNew(pt23-quarter:pt23+quarter) = ...
					( errNew(peakX(pk)-quarter:peakX(pk)+quarter) + ...
					2*errNew(peakX(pk+1)-quarter:peakX(pk+1)+quarter) ) ...
					/3;
			end			
		end
	end
	
	if(voiced == 0)
		disp(['Frame ' num2str(iterator) ' is unvoiced.']);
	end
	
	if(voiced == 0 && m2f == 0 && poleFactor ~= 1) 
		rmsN = norm(err)/sqrt(length(err));
		style1 = rmsN * noiseGain * tukeywin(length(errNew), tukeyshape) .* randn(length(errNew),1);
		errNew = rmsN * noiseGain * tukeywin(length(errNew), tukeyshape) .* randn(length(errNew),1);
	else
		%disp('voiced.');
		style1 = zeros(length(errNew),1); 
		style1(round(peakX * pitchDen/pitchNum)) = errNew(round(peakX * pitchDen/pitchNum));
	end
	
	% Pole scaling 
	a(isnan(a)) = 0;
	[z,porig,k] = tf2zp(1,a);
	p = porig;
	pnew = [];
	if(voiced == 0 || poleFactor == 1) 
		pnew = p;
	else
		if(poleFactor < 1)
			for ii = length(p):-1:1
				if(abs(angle(p(ii))) > .75*pi)
					pnew = [pnew; p(ii)];
					p(ii) = [];
				end
			end
		end
		pnew = [pnew; abs(p) .* exp(1j*angle(p)*poleFactor) .* (abs(angle(p)*poleFactor)<pi)];
		if(poleFactor < 1)
			pnew = [pnew; abs((p(angle(p)==pi))) .* exp(-1j*pi*poleFactor)];
		end
		pnew = pnew(pnew~=0);
	end;
	[num, den] = zp2tf(z, pnew, k);

	if(stylized == 0)
		sh = filter(1, den, errNew);
	else
		sh = filter(1, den, style1);
	end

	errOutput(i:i+length(errNew)-1) = errOutput(i:i+length(errNew)-1) + errNew;
	output(i:i+length(sh)-1) = output(i:i+length(sh)-1) + sh;

	% PLOTS:
	%{
	For male and female:
		-time waveforms
		-LP spectra
		-error signal before and after conversion
		-pole plots
		-spectrograms of original and synthesized signals
	%}
	if (iterator == haltOnIterator) % PLOTTING FINAL RESULTS
	%if(max(sh) > 100) % DEBUGGING
		figure(1);
		subplot(3,1,1); plot(xx); title('Time Waveform'); % TIME WAVEFORM
		%subplot(4,1,2); plot(axx); % AUTOCORRELATION
		subplot(3,1,2); plot(err); title('Error Signal (Original)'); % ERROR SIGNAL
		%hold on; plot(err, 'r'); hold off; % SCALED ERROR SIGNAL
		hold on; plot(peakX,0,'ro'); hold off; % PEAKS
		subplot(3,1,3); plot(sh); title('Error Signal (Scaled and Filtered)');% FILTERED ERROR SIGNAL
		figure(2);
		subplot(2,2,1);
		zplane(z,porig); title('LP Poles (Original)'); % POLE ZERO PLOT (ORIGINAL)
		subplot(2,2,2);
		zplane(z,pnew); title('LP Poles (Scaled)'); % POLE ZERO PLOT (SCALED)
		subplot(2,1,2);
		[h,f] = freqz(1,a,512,fs);
		plot(f,20*log10(abs(h))); xlim([0 fs/2]); title('Linear Prediction Spectral Envelope'); % LP SPECTRA
		[h,f] = freqz(num,den,512,fs);
		hold on; plot(f,20*log10(abs(h)),'r'); legend('Original','Scaled'); xlim([0 fs/2]); hold off; % LP SPECTRA
		if(stylized ~= 0) % STYLIZED SIGNAL
			figure(3);
			plot(style1); title(['Stylized Error Signal (' num2str(stylized) ' peaks per period)']);
		end
		if(halt == 1)
			break;
		end
	end
end

% [output, Noise] = pad(output, Noise);
% output = output + Noise;
output = 0.717 * output / max(abs(output));
output = output(1:length(wav));

figure(4); stairs(mu); axis tight; ylim([-0.05 1.05]); 
hold on; plot(voicedThresh*(ones(length(mu),1)), 'r:'); legend('Detected level','Threshhold','Location','SouthEast'); hold off;
title('Voiced Detection (Entire Signal)'); ylabel('Tonality'); xlabel('Frame'); % VOICE/STOCHASTIC DETECTOR
% plot(output);
figure(5); N=512;
subplot(2,1,1); specgram(wav,N,fs,hamming(N/4),round(0.9*N/4));xlabel('time, s');
title('Spectrogram (Original)');
subplot(2,1,2); specgram(output,N,fs,hamming(N/4),round(0.9*N/4));xlabel('time, s');
title('Spectrogram (Synthesized)');

%findfigs;
soundsc(output, fs);
if(stylized == 0) sty = []; else  sty = ['_stylized' num2str(stylized)]; end
if(child == 1) ch = '_child'; else ch = '_female'; end
if(m2f == 0) filetitle = [filename(1:end-4) sty '_to_male'];
else filetitle = [filename(1:end-4) sty '_to' ch]; end

set(figure(1), 'PaperPosition', [0 0 6 4]);
set(figure(2), 'PaperPosition', [0 0 6 4]);
set(figure(3), 'PaperPosition', [0 0 6 4]);
set(figure(4), 'PaperPosition', [0 0 6 4]);
set(figure(5), 'PaperPosition', [0 0 6 4]);
print(figure(1),[filetitle '_analysis.png'],'-dpng','-r100');
print(figure(2),[filetitle '_lp.png'],'-dpng','-r100');
if(stylized ~= 0) print(figure(3), [filetitle '_stylized.png'],'-dpng','-r100'); end
print(figure(4), [filetitle '_voiced.png'],'-dpng','-r100');
print(figure(5), [filetitle '_specgram.png'],'-dpng','-r100');
wavwrite(output, fs, [filetitle '.wav']);
if(poleFactor == 1) 
	errOutput = .717 * errOutput / max(abs(errOutput));
	wavwrite(errOutput, fs, [filetitle '_error.wav']);
end
return;

%% 
% Caruso restoration

clear all; close all; clc;
[caruso, fs] = wavread('caruso_short.wav');
[pavarotti, fs] = wavread('pavarotti.wav');

L = 1024;
alsCaruso = averageLogSpectrum(caruso,L);
subplot(2,1,1);
plot(abs(alsCaruso(L/2:end)));
alsPavarotti = averageLogSpectrum(pavarotti,L);
hold on; plot(abs(alsPavarotti(L/2:end)),'r'); hold off;
subplot(2,1,2);
alsDiff = abs(alsCaruso)-abs(alsPavarotti); % (L/2:end)
plot(alsDiff,'m');
asDiff = exp(alsDiff);
figure(2);
subplot(3,1,1); plot(asDiff);
H = 1./asDiff;
subplot(3,1,2); plot(H);
h = real(ifft(H));
subplot(3,1,3); plot(h);

h = h(1:L/2);
carusoNew = conv(caruso,h);

	% magnitude spectrum plots
	figure(3);
	subplot(2,1,1); fftplot(caruso,fs);
	hold on; fftplot(carusoNew,fs,'r'); hold off;
	title('Caruso Magnitude Spectrum');
	legend('Original','Restored');
	subplot(2,1,2); fftplot(pavarotti,fs);
	title('Pavarotti Magnitude Spectrum');


soundsc([carusoNew(fs*2:fs*4); caruso(fs*4:fs*6); carusoNew(fs*6:fs*8); caruso(fs*8:fs*10) ],fs);

