% FFTconv.m
% DSP - eMe, 2010 - martin rocamora
% ejemplo de calculo de convolucion usando la FFT
close all
clear all

% entrada de audio (mono)
[input,fs,nbits] = wavread('ohwhere_22050.wav');

% respuesta al impulso de sala
[h fs_h] = wavread('FOLIP-12.wav');
% izquierda y derecha
h_l = h(:,1);
h_r = h(:,2);

% largo de cada senial
largo_input = length(input);  
largo_h     = length(h_l);

% puntos de ventana de analisis
N = 2048;

% puntos de la FFT, se agregan suficientes ceros 
% para evitar convolucion circular
N_fft = N + largo_h - 1; 
 
% transformada de las respuestas al impulso
Fh_l = fft(h_l,N_fft);    
Fh_r = fft(h_r,N_fft);

% numero de tramas
num_ventantas = floor(largo_input/N)

% senial de salida
output_l = zeros(num_ventantas * N + largo_h-1, 1);
output_r = zeros(num_ventantas * N + largo_h-1, 1);

for i = 1:num_ventantas
   
   % puntos de inicio y fin de la trama 
   ini = (i-1) * N + 1;
   fin = ini + N - 1;
    
   % trama de senial 
   trama = input(ini:fin);
   
   % Fourier de la trama
   Ftrama = fft(trama, N_fft);
   
   % producto de espectros
   F_l = Ftrama .* Fh_l; 
   F_r = Ftrama .* Fh_r;
   
   % se antitransforma y se toma parte real (parte compleja deberia ser cero)
   trama_l = real(ifft(F_l));	
   trama_r = real(ifft(F_r));
   
   % se almacena y solapa
   output_l(ini:(ini+N_fft-1),1) = output_l(ini:(ini+N_fft-1),1) + trama_l;
   output_r(ini:(ini+N_fft-1),1) = output_r(ini:(ini+N_fft-1),1) + trama_r;
      
   %i
      
end

% normalizar amplitud (a la misma que la entrada)
max_amp = max(max(abs(output_l)), max(abs(output_r)));
max_input = max(abs(input));
output_l = output_l * max_input / max_amp;
output_r = output_r * max_input / max_amp;

% guardar salida
output = [output_l output_r];          
wavwrite(output, fs, nbits,'FFTconv');


%% %%%%%%%% GRAFICAS %%%%%%%%%

t = (0:(num_ventantas * N + largo_h-1)-1)/fs;

figure
subplot(311)
plot(t(1:largo_input), input, 'k')
ylabel('amplitud')
title('Senial de entrada y respuestas al impulso')
subplot(312)
plot(t(1:largo_h), h_l, 'k')
ylabel('amplitud')
subplot(313)
plot(t(1:largo_h), h_r, 'k')
ylabel('amplitud')
xlabel('tiempo (s)')

figure
subplot(211)
plot(t, output_l, 'k')
ylabel('amplitud')
title('Seniales de salida')
subplot(212)
plot(t, output_r, 'k')
ylabel('amplitud')
xlabel('tiempo (s)')



