clc;clear all;close all
load resultados_sin_cabeza
%Cargamos el programa de cálculo del tiempo de reverberación por tercios de
%octava
fs=44100;
for ii=1:16
[rt0(:,ii),output2,yy0(:,:,ii)]=rt(salida0(:,ii),44100,1,-5,-15,1);axis( [1  2.5  -60 0]);  
[rt1(:,ii),output2,yy1(:,:,ii)]=rt(salida1(:,ii),44100,1,-5,-15,1);axis( [1  2.5  -60 0]);
[rt2(:,ii),output2,yy2(:,:,ii)]=rt(salida2(:,ii),44100,1,-5,-15,1);axis( [1  2.5  -60 0]);
end
%Se realiza un for para promediar las caidas de los 16 micrófonos colocados
for ii=1:18
ayy0(:,:)=yy0(:,ii,:);
byy0(:,ii)=mean(ayy0');
end
for ii=1:18
ayy1(:,:)=yy1(:,ii,:);
byy1(:,ii)=mean(ayy1');
end
for ii=1:18
ayy2(:,:)=yy2(:,ii,:);
byy2(:,ii)=mean(ayy2');
end
%PLOTEA LAS CAIDAS DE CADA SALA
%SALA 0
plotparameter=1;
bands=1;
if plotparameter==1
   plot((0:length(byy0')-1)/fs,byy0')
   
   if bands==1
      legend('100','125','160','200','250','315','400','500','630','800','1000','1250','1600','2000','2500','3150','4000','5000');
   else
      legend('125','250','500','1000','2000','4000');
   end
axis( [1  2.1  -60 0]) 
title('Curvas de caida del nivel de presión sonora. Sala sin difusores')
xlabel('Tiempo')
ylabel('dB')
end

%SALA 0
figure
plotparameter=1;
bands=1;
if plotparameter==1
   plot((0:length(byy1')-1)/fs,byy1')
   
   if bands==1
      legend('100','125','160','200','250','315','400','500','630','800','1000','1250','1600','2000','2500','3150','4000','5000');
   else
      legend('125','250','500','1000','2000','4000');
   end
axis( [1  2.1  -60 0]) 
title('Curvas de caida del nivel de presión sonora. Sala con cristales de sonido convencionales')
xlabel('Tiempo')
ylabel('dB')
end

%SALA 0
figure
plotparameter=1;
bands=1;
if plotparameter==1
   plot((0:length(byy2')-1)/fs,byy2')
   
   if bands==1
      legend('100','125','160','200','250','315','400','500','630','800','1000','1250','1600','2000','2500','3150','4000','5000');
   else
      legend('125','250','500','1000','2000','4000');
   end
axis( [1  2.1  -60 0]) 
title('Curvas de caida del nivel de presión sonora. Sala con cristales de sonido optimizados')
xlabel('Tiempo')
ylabel('dB')
end
%Ajuste de rectas

%sala 0
A0=((0:length(byy0')-1)/fs)';
A0=A0(50000:100000);
B0=byy0(50000:100000,:);
for ii=1:18
  C0(ii,:)=polyfit(A0',B0(:,ii)',1);
  y0(:,ii)=C0(ii,1).*A0+C0(ii,2);
  coef0(ii)=max(xcorr(y0(:,ii),B0(:,ii),'coeff'));
  cor0(ii)=max(corr(y0(:,ii),B0(:,ii)));
end

%sala 1
A1=((0:length(byy1')-1)/fs)';
A1=A1(50000:76000);
B1=byy1(50000:76000,:);
for ii=1:18
  C1(ii,:)=polyfit(A1',B1(:,ii)',1);
  y1(:,ii)=C1(ii,1).*A1+C1(ii,2);
  coef1(ii)=max(xcorr(y1(:,ii),B1(:,ii),'coeff'));
  cor1(ii)=max(corr(y1(:,ii),B1(:,ii)));
end

%sala 2
A2=((0:length(byy2')-1)/fs)';
A2=A2(50000:76000);
B2=byy2(50000:76000,:);
for ii=1:18
  C2(ii,:)=polyfit(A2',B2(:,ii)',1);
  y2(:,ii)=C2(ii,1).*A2+C2(ii,2);
  coef2(ii)=max(xcorr(y2(:,ii),B2(:,ii),'coeff'));
  cor2(ii)=max(corr(y2(:,ii),B2(:,ii)));
end
%ploteando ejemplo de ajuste de rectas
subplot(3,1,1)
plot(A0,y0(:,4),'g',A0,B0(:,4),'r','LineWidth',2)

title('Ajuste por mínimos cuadrados en sala sin difusores')
xlabel('Tiempo')
ylabel('dB')
legend('Recta','Curva caida 200Hz')

subplot(3,1,2)
plot(A1,y1(:,4),'b',A1,B1(:,4),'p','LineWidth',2)

title('Ajuste por mínimos cuadrados en sala sin difusores')
xlabel('Tiempo')
ylabel('dB')
legend('Recta','Curva caida 200Hz')

subplot(3,1,3)
plot(A2,y2(:,4),'k',A2,B2(:,4),'r','LineWidth',2)

title('Ajuste por mínimos cuadrados en sala sin difusores')
xlabel('Tiempo')
ylabel('dB')
legend('Recta','Curva caida 200Hz')

%escritura de resultados en EXCEL
xlswrite('Resultados.xlsx',rt0,'EDT');
xlswrite('Resultados.xlsx',rt1,'EDT','A20');
xlswrite('Resultados.xlsx',rt2,'EDT','A40');

% xlswrite('Resultados.xlsx',coef0','Curvas');
% xlswrite('Resultados.xlsx',cor0','Curvas','B1');
% xlswrite('Resultados.xlsx',coef1','Curvas','C1');
% xlswrite('Resultados.xlsx',cor1','Curvas','D1');
% xlswrite('Resultados.xlsx',coef2','Curvas','E1');
% xlswrite('Resultados.xlsx',cor2','Curvas','F1');


return
A=((0:length(byy0')-1)/fs)';
B=ayy0(:,1);
B1=B(50000:100000);
A1=A(50000:100000);
C=polyfit(A1',B1',1);

y=C(1).*A1+C(2);

plot(A1,y);hold on;
plot(A1,B1);

max(xcorr(y,B1,'coeff'));
