% test_pv - test phase vocoder
% 
close all
clear all

% archivo de audio
%[d,fs] = wavread('ohwhere_22050.wav');
[d,fs] = wavread('love.wav');

%% ESTIRAMIENTO TEMPORAL

% largo ventana
N = 1024;
% factor de estiramiento temporal
ts = 0.75;

% phase vocoder
y = pvoc(d, ts, N);

wavwrite(y,fs,'salida_stretch.wav')

%% ARMONIZACION: CORRIMIENTO DE ALTURA Y SUMA

% factor de estiramiento temporal
ta = 0.8; % 0.8 = 4/5
p  =   4;
q  =   5;

% ta = 1.6; % 1.6 = 8/5
% p  =   8;
% q  =   5;

% ta = 0.5; % 0.5 = 1/2
% p  =   1;
% q  =   2;

% phase vocoder
e = pvoc(d, ta, N);
% remuestreo
f = resample(e, p, q); 

% suma de ambas señales
g = d(1:length(f)) + f;

wavwrite(g,fs,'salida_harmonization.wav')

%% GRAFICAS

[Sd, fd, td] = specgram(d, N, fs, hanning(N), N-N/4);
[Sy, fy, ty] = specgram(y, N, fs, hanning(N), N-N/4);
[Sf, ff, tf] = specgram(f, N, fs, hanning(N), N-N/4);

fmax = 2000;
indd_max = find(fd > fmax, 1, 'first') - 1;
indy_max = find(fy > fmax, 1, 'first') - 1;
indf_max = find(ff > fmax, 1, 'first') - 1;

ttd = ((1:length(d))-1)/fs;
tendd = find(ttd > td(end), 1, 'first') - 1 + N/2; 
tty = ((1:length(y))-1)/fs;
tendy = find(tty > ty(end), 1, 'first') - 1 + N/2; 


ylims = [-0.3 0.3];

figure;
subplot(211)
plot(ttd(1:tendd), d(1:tendd), 'k')
title('Original')
xlabel('Tiempo (s)')
ylabel('Amplitud')
%xlim([0 ttd(tendd)])
xlim([0 tty(tendy)])
ylim(ylims)
subplot(212)
imagesc(td,fd(1:indf_max),20*log10(abs(Sd(1:indd_max,:))))
colormap(abs(1-gray))
axis xy
xlim([0 tty(tendy)])
title('Espectrograma')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

print -depsc 'pv_original.eps'

figure;
subplot(211)
plot(tty(1:tendy), y(1:tendy), 'k')
title('Escalado temporal')
xlabel('Tiempo (s)')
ylabel('Amplitud')
xlim([0 tty(tendy)])
ylim(ylims)
subplot(212)
imagesc(ty,fy(1:indf_max),20*log10(abs(Sy(1:indd_max,:))))
colormap(abs(1-gray))
axis xy
title('Espectrograma')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

print -depsc 'pv_stretching.eps'

figure;
subplot(211)
imagesc(td,fd(1:indd_max),20*log10(abs(Sd(1:indd_max,:))))
colormap(abs(1-gray))
axis xy
title('Espectrograma original')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')
subplot(212)
imagesc(tf,ff(1:indf_max),20*log10(abs(Sf(1:indf_max,:))))
colormap(abs(1-gray))
axis xy
title('Espectrograma de trasposicion')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

print -depsc 'pv_trasposicion.eps'
