%% Phase Vocoder
% Trabalho 01 de Processameto Digital de Sinais 2
% Gisele Tam - GRR20102167
% 
% Este código foi baseado no código escrito por D.P.W. Ellis em 2002, que
% está disponível em http://labrosa.ee.columbia.edu/matlab/pvoc/

[original,fs] = wavread('singing.wav');
sound(original,fs)
display('Som original')

windowLen = 1024; %tamanho da janela; 
nJump = (windowLen)/4; %passo da janela

% Parâmetros para compressão/dilatamento no tempo e mudança de frequência
% Para mais agudo: high>low; 
% Para mais grave: high<low.
low = 3;
high = 4;
%% Aplicação da STFT (Short-time Fourier Transform)

sizeO = length(original);
vAnalise = original(:,1)';

mSTFT = zeros((1+windowLen/2),1+floor(sizeO-windowLen)/nJump);
win = (hann(windowLen))'; %janela Hanning
col = 1; %inicialização do valor da coluna
for ind = 0:nJump:(sizeO-windowLen)
    aux = win.*vAnalise((ind+1):(ind+windowLen));
    aux = fft(aux);
    mSTFT(:,col) = (aux(1:(1+windowLen/2)))'; %usa metade para retirar os pontos redundantes
    col = col+1;
end


%% Novo spectograma

[r,c] = size(mSTFT);
factor = low/high;
vChange = 0:factor:(c-2);
N = 2*(r-1);

mNew = zeros(r,length(vChange)); %inicialização do vetor de interpolação

% Acumulador de fase
ph = angle(mSTFT(:,1));

% Matriz auxiliar equivalente a mNew com uma coluna de segurança de zeros
% no final caso seja um número ímpar de colunas
mAux = [mSTFT, zeros(r,1)];

col = 1; %auxiliar para coluna no loop
for fac = vChange
    mAuxC = mAux(:,floor(fac)+[1 2]);
    ff = fac - floor(fac);
    mMag = (1-ff)*abs(mAuxC(:,1)) + ff*(abs(mAuxC(:,2)));
    % Cálulo do avanço de fase
    dp = angle(mAuxC(:,2)) - angle(mAuxC(:,1));
    % Deixar os valores entre -pi e pi
    dp = dp -2*pi*round(dp/(2*pi));
    % Nova matriz
    mNew(:,col) = mMag.*exp(j*ph);
    % Acumulação da fase 
    ph = ph + dp;
    col = col+1;
end

%% Aplicação da ISTFT (STFT Inversa)

[rows, cols] = size(mNew);

winISTF = win;
lengthISTF = windowLen + (cols-1)*nJump;
vISTF = zeros(1,lengthISTF);

for ind = 0:nJump:(nJump*(cols-1))
    i_aux = (mNew(:,1+ind/nJump))';
    i_aux = [i_aux, conj(i_aux([(windowLen/2):-1:2]))];
    vecI = real(ifft(i_aux));
    vISTF((ind+1):(ind+windowLen)) = vISTF((ind+1):(ind+windowLen))+vecI.*winISTF;
end

vISTF = vISTF;

%% 
% Comprimido/dilatado no tempo
sound(vISTF,fs)
display('Som após PhaseVocoder')

% Mudança de frequêmcia
pShift = resample(vISTF,low,high);
sound(pShift,fs)
display('Som após mudança de frequência')

