%% Project 2

clear all;
close all;

%% Part 1 & 2

% Use a modern recording of Pavarotti provided on Bb (or similar) to estimate the average,
% long-term spectrum of the tenor?s voice. Plot the log-magnitude spectrum. Plot the
% principal phase spectrum.
% 
% Develop a modulo-2? phase unwrapping method and use it to produce the actual phase
% spectrum of the modern tenor?s voice.

N = 1024;                                               % Window Size and FFT length

[pavarotti, fs] = audioread('PavarottiShort.wav');      % Read the Pavarotti audio file
pavarotti = pavarotti';                                 % Invert the vector
pavarotti = pavarotti./max(abs(pavarotti));             % Normalize the audio
pavFFT = averageSpectrum(pavarotti, N, N);              % Get it's average spectrum

pavMagDB = mag2db(abs(pavFFT));                         % Log Magnitude Spectrum
pavPhase = angle(pavFFT);                               % Principle Phase Spectrum
pavPhaseUnwrapped = unwrap(pavPhase);                   % Unwrapped Phase Spectrum

pavIFFT = ifft(pavFFT, N);                              % Take the IFFT for LPA
pavLPA = linearPredictionAnalysis(pavIFFT, fs, N);      % Do LPA

freq = 0:fs/N:fs/2-fs/N;                                % Frequency Vector

fig = figure(1);
set(fig, 'Position', [100 100 1000 800]);

subplot(3,1,1)
plot(freq, pavMagDB(1:N/2));                            % Plot the Log Mag Spectrum
hold on;
plot(freq, pavLPA(1:N/2));                              % and superimpose the LPA curve
title('Pavarotti Log-Magnitude Spectrum')
axis([freq(1) freq(end) -128 6]);
xlabel('Frequency');
ylabel('Magnitude_{dB}');

subplot(3,1,2)
plot(freq, pavPhase(1:N/2));                            % Plot the Principle Phase Spectrum
title('Pavarotti Wrapped Phase Spectrum')
axis([freq(1) freq(end) -4 4]);
xlabel('Frequency');
ylabel('Phase Shift_{radians}');

subplot(3,1,3)
plot(freq, pavPhaseUnwrapped(1:N/2));                   % Plot the Unwrapped Phase Spectrum
title('Pavarotti Unwrapped Phase Spectrum')
axis([freq(1) freq(end) -40 40]);
xlabel('Frequency');
ylabel('Phase Shift_{radians}');

saveas(fig, 'fig1_Project2', 'png');


%% Part 3

% Estimate the average, long-term spectrum of the available Caruso recording. Plot the logmagnitude,
% the principal phase, and the unwrapped phase spectrum.

[caruso, fs2] = audioread('CarusoShort.wav');           % Read the Caruso audio file
caruso = caruso';                                       % Invert the vector
carusoL = caruso(1,:);                                  % Use only the left channel
carusoL = carusoL./max(abs(carusoL));                   % Normalize the audio
carFFT = averageSpectrum(carusoL, N, N);                % Get it's average spectrum

carMagDB = mag2db(abs(carFFT));                         % Log Magnitude Spectrum
carPhase = angle(carFFT);                               % Principle Phase Spectrum
carPhaseUnwrapped = unwrap(carPhase);                   % Unwrapped Phase Spectrum

carIFFT = ifft(carFFT, N);                              % Take the IFFT for LPA
carLPA = linearPredictionAnalysis(carIFFT, fs, N);      % Do LPA

freq = 0:fs/N:fs/2-fs/N;                                % Frequency Vector

fig = figure(2);
set(fig, 'Position', [100 100 900 700]);

subplot(3,1,1)
plot(freq, carMagDB(1:N/2));
hold on;
plot(freq, carLPA(1:N/2));
title('Caruso Log-Magnitude Spectrum');
axis([freq(1) freq(end) -128 6]);
xlabel('Frequency');
ylabel('Magnitude_{dB}');

subplot(3,1,2)
plot(freq, carPhase(1:N/2));
title('Caruso Wrapped Phase Spectrum');
axis([freq(1) freq(end) -4 4]);
xlabel('Frequency');
ylabel('Phase Shift_{radians}');

subplot(3,1,3)
plot(freq, carPhaseUnwrapped(1:N/2));
title('Pavarotti Unwrapped Phase Spectrum');
axis([freq(1) freq(end) -160 10]);
xlabel('Frequency');
ylabel('Phase Shift_{radians}');

saveas(fig, 'fig2_Project2', 'png');

%% Part 4 

carIFFT = ifft(carFFT, N);                              % Take the IFFT of the two
pavIFFT = ifft(pavFFT, N);                              % FFTs (Caruso & Pavarotti)

carCepstral = cceps(carIFFT);                           % Go to cepstral domain
pavCepstral = cceps(pavIFFT);

carCepstral = cat(2, carCepstral(N/2+1:N), carCepstral(1:N/2));     % Rearrange cepstrums
pavCepstral = cat(2, pavCepstral(N/2+1:N), pavCepstral(1:N/2));

transCepstral = carCepstral - pavCepstral;              % Subtract Pav from Car to get the
                                                        % cepstrum of the transducer

n = 0:1:N-1;                                            % Quefrency vector

fig = figure(3);
set(fig, 'Position', [100 100 1000 800]);

subplot(3,1,1)
plot(n, carCepstral);                                   % Plot all of the cepstrums
title('Caruso Complex Ceptrum');
axis([n(1) n(end) min(carCepstral) max(carCepstral)]);
xlabel('Quefrency');
ylabel('Amplitude');

subplot(3,1,2)
plot(n, pavCepstral);
title('Pavarotti Complex Ceptrum');
axis([n(1) n(end) min(carCepstral) max(carCepstral)]);
xlabel('Quefrency');
ylabel('Amplitude');

subplot(3,1,3)
plot(n, transCepstral);
title('Transducer Complex Ceptrum (Caruso - Pavarotti)');
axis([n(1) n(end) min(carCepstral) max(carCepstral)]);
xlabel('Quefrency');
ylabel('Amplitude');

saveas(fig, 'fig3_Project2', 'png');

transducerIR = icceps(transCepstral);                   % Get the Impulse Response of the transducer...

windowedTransducerIR = transducerIR .* hamming(length(transducerIR))';  % Window the transducerIR

transducerFFT = fft(windowedTransducerIR, N);           % Take the FFT of the transducerIR

transMagDB = mag2db(abs(transducerFFT));                % Log Magnitude Spectrum
transPhase = angle(transducerFFT);                      % Wrapped Phase Spectrum
transPhaseUnwrapped = unwrap(transPhase);               % Unwrapped Phase Spectrum

[transLPA, transLPAcoeffs, transLPAorder] = linearPredictionAnalysis(windowedTransducerIR, fs, N); % Do LPA

fig = figure(4);
set(fig, 'Position', [100 100 1000 800]);

subplot(2,2,1)
plot(n, windowedTransducerIR);
title('Windowed Transducer Impulse Response');
axis([n(1) n(end) min(windowedTransducerIR) max(transducerIR)]);
xlabel('Quefrency');
ylabel('Amplitude');

subplot(2,2,2);
plot(freq, transMagDB(1:N/2));
hold on;
plot(freq, transLPA(1:N/2));
title('Windowed Transducer Log-Magnitude Response');
axis([freq(1) freq(end) -32 64]);
xlabel('Quefrency');
ylabel('Amplitude');

subplot(2,2,3)
plot(freq, transPhase(1:N/2));
title('Windowed Transducer Wrapped Phase Response');
axis([n(1) n(end) -4 4]);
xlabel('Quefrency');
ylabel('Amplitude');

subplot(2,2,4)
plot(freq, transPhaseUnwrapped(1:N/2));
title('Windowed Transducer Unwrapped Phase Response');
axis([n(1) n(end) -20 20]);
xlabel('Quefrency');
ylabel('Amplitude');

saveas(fig, 'fig4_Project2', 'png');

%% Part 5

fig = figure(5);
set(fig, 'Position', [100 100 1000 800]);
poles = roots(transLPAcoeffs);
zplane([], poles);
title('Z Plane representation of the Transducer"s All-Pole System Function');


saveas(fig, 'fig5_Project2', 'png');