clear all;close all
% MÃ©todo de deteccion de daÃ±o basado en transformada wavelet desarrollado
% por M. Algaba, M. SolÃ­s y P. GalvÃ­n (Grupo Estructuras ETSI Univ.
% Sevilla)
%
%
% VERSIÃ“N V10 INCLUYE CONSIDERACIÃ“N DE COORDENADAS DE PUNTOS DE MEDIDA Y
% TOMANDO SIEMPRE MUESTRAS EXPERIMENTALES EN SECCIONES EXTREMAS DE LA VIGA
% (AFECTA A INTERPOLACIONES, EXTENSIONE)
% TAMBIEN CONSIDERA LOS 65 PUNTOS EN VEZ DE 64 (CON CWT NO AFECTA EL TAMAÃ‘O
% DIADICO, AL MENOS CON DAUBECHIES 2)
%
% DefiniciÃ³n de nÃºmeros de muestras para reducciÃ³n muestral (si es igual
% que el tamaÃ±o del vector original, 65, entonces no se aplicarÃ¡ despuÃ©s la
% reducciÃ³n)

%tam_reduc=15;%128,64,32;
% DefiniciÃ³n de numero demuestras para la interpolaciÃ³n (este serÃ¡ el
% nÃºmero de muestras comprendidas en la longitud de la viga)
tam_reduc=128;
tam_interpol=65;%128,65,64 
%
% DefiniciÃ³n del vector de posiciÃ³n de los acelerÃ³metros
%
for i=1:128 %65,128
    x(i)=(i-1)*.0101;%0.02,0.0101
end
%x=cat(2,x(1:52),x(54:65));
%x=cat(2,x(1:32),x(34:65));

opcion=menu('Cargar archivos con las propiedades dinámicas de la barra SIN daño','Buscar archivo con las frecuencia y amortiguamiento',...
 'Salir');

%Cargar archivo con los vectores de frecuencia y amortiguación
if opcion==1
uiload
else %Salida del programa
break
end


opcion=menu('Cargar archivos con las propiedades dinámicas de la barra SIN daño','Buscar archivo con los modos reales',...
 'Salir');

%Cargar archivo con la matriz de modos
if opcion==1
uiload
else %Salida del programa
break
end
matriz_modo_intacto_exp=mod_real;
frec_intacta=freq_sel;
%matriz_modo_intacto_exp=cat(1,mod_real(1:52,:),mod_real(54:65,:));
clear mod_real;
% % 
opcion=menu('Cargar archivos con las propiedades dinámicas de la barra CON daño','Buscar archivo con las frecuencia y amortiguamiento',...
 'Salir');

%Cargar archivo con los vectores de frecuencia y amortiguación
if opcion==1
uiload
else %Salida del programa
break
end


opcion=menu('Cargar archivos con las propiedades dinámicas de la barra CON daño','Buscar archivo con los modos reales',...
 'Salir');

%Cargar archivo con la matriz de modos
if opcion==1
uiload
else %Salida del programa
break
end


matriz_modo_defecto_exp=mod_real;
frec_defecto=freq_sel;
% matriz_modo_defecto_exp=mod_real(1:65,:);
% matriz_modo_defecto_exp=cat(1,mod_real(1:52,:),mod_real(54:65,:)); %Se
% elimina picotazo de muestra 53
%
clear modes_real;

%
% Cargamos los modos de la viga que corresponda y se pone en la matriz de
% modos intacta y daÃ±ada  (se toman 64 puntos para
% mantener la seÃ±al diadica)
%
%load('Modes_CSI/Intacta_conf_ALL_proc_dec2_vt8_csi_i16_ord60_modes_real.mat');
%matriz_modo_intacto_exp=mod_real(1:65,:);
%frec_intacta=[415.65
%1032.70
%1786.75
%2581.50
%3366.63];
%clear mod_real;
% % 
% load('Modes_CSI/05Lmax_conf_all_CSI_i16_ord60_modes_real.mat');
% matriz_modo_defecto_exp=modes_real(1:65,:);
% frec_defecto=[300.96
% 1027.35
% 1473.96
% 2563.38
% 3195.89];
% clear modes_real;
% 
% load('Modes_CSI/05L_med_conf_ALL_proc_dec2_vt8_csi_i16_ord60_modes_real.mat');
% matriz_modo_defecto_exp=mod_real(1:65,:);
% frec_defecto=[362.28
% 1030.15
% 1634.17
% 2576.63
% 3253.80];
% clear mod_real;
% % % 
% load('Modes_CSI/05Lleve_ALL_but_conf1_ev_2-3_proc_dec2_vt32_csi_i16_ord60_modes_real.mat');
% matriz_modo_defecto_exp=mod_real(1:65,:);
% matriz_modo_defecto_exp=cat(1,mod_real(1:52,:),mod_real(54:65,:)); %Se elimina picotazo de muestra 53
% clear mod_real
% clear matriz_modo_intacto_exp
% load('Modes_CSI/Intacta_conf_ALL_proc_dec2_vt8_csi_i16_ord60_modes_real.mat');
% matriz_modo_intacto_exp=cat(1,mod_real(1:51,:),mod_real(54:65,:));
% x=cat(2,x(1:51),x(54:65));
% if tam_reduc==65 % Al no poder haber 65 muestras se toman 64 y habrÃ¡ que hacer reduccion muestral forzosamente
%     tam_reduc=63
% end
% frec_defecto=[407.57
% 1033.32
% 1781.58
% 2566.85
% 3354.80];
% clear mod_real
% % % 
% load('Modes_CSI/025Lmax_conf_All_but_conf3-4-5-5-7-ev1-3-1-2-1_proc_dec2_vt32_csi_i16_ord60_modes_real.mat');
% matriz_modo_defecto_exp=mod_real(1:65,:);
% frec_defecto=[364.75
% 822.99
% 1557.21
% 2528.65
% 3292.95];
% clear mod_real;
% % 
% load('Modes_CSI/025Lmed_conf_ALL_proc_dec2_vt32_csi_i16_ord60_modes_real.mat');
% matriz_modo_defecto_exp=mod_real(1:65,:);
% frec_defecto=[397.27
% 932.71
% 1663.19
% 2556.50
% 3354.27];
% clear mod_real;
% % 
% load('Modes_CSI/025L_leve_conf_ALL_but_conf_4-4-5_ev_2-4-5_proc_dec2_vt32_csi_i16_ord60_modes_real.mat');
% mod_real(:,6)=[];
% mod_real(:,4)=[];% se borran 4Âª y 6Âª columna porque modos 4 y 5 estÃ¡n repetidos en el fichero original. Esto no afecta aresultados si solo se trabaja con 3 modos
% matriz_modo_defecto_exp=cat(1,mod_real(1:32,:),mod_real(34:65,:)); %elimina muestra 33, que tiene picotazo de ruido, y hace lo mismo con el vector intacto
% clear mod_real;
% clear matriz_modo_intacto_exp
% load('Modes_CSI/Intacta_conf_ALL_proc_dec2_vt8_csi_i16_ord60_modes_real.mat');
% matriz_modo_intacto_exp=cat(1,mod_real(1:32,:),mod_real(34:65,:));
% clear mod_real;
% x=cat(2,x(1:32),x(34:65));
% if tam_reduc==65 % Al no poder haber 65 muestras se toman 64 y habrÃ¡ que hacer reduccion muestral forzosamente
%     tam_reduc==64
% end
% frec_defecto=[415.01
% 1020.23
% 1772.57
% 2580.6
% 3371.16];

%load('Modes_CSI/borde_conf_all_event_all_but conf3_ev134_proc_dec2_vt32_csi_i16_ord60_modes_real.mat');
%matriz_modo_defecto_exp=mod_real(1:65,:);
%frec_defecto=[417.80
%1034.59
%1790.30
%2564.90
%3096.85];
%clear mod_real;
% 

% 
% Para reducir el tamaÃ±o muestral.

% 
tam=size(matriz_modo_intacto_exp);
tam_modo=tam(1);
num_modos=tam(2);
% 
% 
if tam_reduc==tam_modo
        matriz_modo_intacto_exp=matriz_modo_intacto_exp;
        matriz_modo_defecto_exp=matriz_modo_defecto_exp;
else
        [matriz_modo_intacto_exp,x_reduc] = reduccion_muestral_v4 (matriz_modo_intacto_exp,x,tam_reduc);
        [matriz_modo_defecto_exp,x_reduc] = reduccion_muestral_v4 (matriz_modo_defecto_exp,x,tam_reduc);
        x=x_reduc;
end


% En primer lugar comprobamos que estï¿½n en FASE los modos. Para ello
% realizamos el producto escalar de los modos (intacto y daï¿½ado) y si se
% obtiene un valor negativo, habra que multiplicar por (-1) y en caso
% contrario se deja como estï¿½. Se decide cambiar de fase el modo daï¿½ado.
%
for i=1:num_modos
    if matriz_modo_intacto_exp(:,i)'* matriz_modo_defecto_exp(:,i)<0
        matriz_modo_defecto_exp(:,i)=-matriz_modo_defecto_exp(:,i);
    end
end
% % Se normalizan los modos a la unidad
% % %
% 
% for i=1:5
%          matriz_modo_intacto_exp(:,i)=matriz_modo_intacto_exp(:,i)/norm(matriz_modo_intacto_exp(:,i));
%         matriz_modo_defecto_exp(:,i)=matriz_modo_defecto_exp(:,i)/norm(matriz_modo_defecto_exp(:,i));
%     
% end
%
% % %
% %  Determinamos el nï¿½mero de modos considerados y el tamaï¿½o
% muestral de los mismos.
%
tam=size(matriz_modo_intacto_exp);
tam_muestral=tam(1);
%
%

%
%Se suaviza mediante interpolaciï¿½n, en el caso de querer introducir muestras
%intermedias entre las muestras reducidas o porque las muestras
%experimentales no estÃ©n equiespaciadas y distribuidas uniformemente a lo
%largo de la viga
% Si el nÃºmero de muestras de la seÃ±al reducida es igual al de la seÃ±al
% original, no se realiza la interpolaciÃ³n y se realiza un suavizado del
% modo mediante curve-fitting

%
for i=1:num_modos
        matriz_modo_intacto_ext(:,i)=wextend(1,'asymw',matriz_modo_intacto_exp(:,i),(tam_muestral-1),'b');
        matriz_modo_defecto_ext(:,i)=wextend(1,'asymw',matriz_modo_defecto_exp(:,i),(tam_muestral-1),'b');
end
    x_ext=wextend(1,'asymw',x,tam_muestral-1,'b');

if tam_reduc==tam_modo;
   [matriz_modo_intacto_soft_ext,matriz_modo_defecto_soft_ext]=curve_fitting_v3(matriz_modo_intacto_ext,...
    matriz_modo_defecto_ext,x_ext);

else 
     
    for i=1:num_modos
          [matriz_modo_intacto_soft_ext(:,i),x_soft_ext]=interpolacion_spline_CWT_v4(matriz_modo_intacto_ext(:,i),x_ext,3*tam_interpol-2);
          [matriz_modo_defecto_soft_ext(:,i),x_soft_ext]=interpolacion_spline_CWT_v4(matriz_modo_defecto_ext(:,i),x_ext,3*tam_interpol-2);
    end
   x_ext=x_soft_ext;
end
%
%

%
% Calculamos las diferencias de los modos ponderadas por sus frecuencias
% naturales. 
%
for j=1:num_modos
    for i=1:length(matriz_modo_defecto_soft_ext(:,1))
        dif_modos(i,j)=(matriz_modo_defecto_soft_ext(i,j)-matriz_modo_intacto_soft_ext(i,j))*...
       (1-frec_intacta(j)/frec_defecto(j))^2;
       
    end
end
%
% Calculamos la transformada wavelet continua para cada uno de las diferencias 
% de los modos introducidos en la matriz
%
N_scale=log10(length(matriz_modo_defecto_soft_ext(:,1)))/log10(2);
coefs_dif_modo1 = CWT(dif_modos(:,1),1:N_scale,'db2');
coefs_dif_modo2 = CWT(dif_modos(:,2),1:N_scale,'db2');
coefs_dif_modo3 = CWT(dif_modos(:,3),1:N_scale,'db2');
coefs_dif_modo4 = CWT(dif_modos(:,4),1:N_scale,'db2');
coefs_dif_modo5 = CWT(dif_modos(:,5),1:N_scale,'db2');
%
% Ploteamos la matriz de diferencias de modos
%
% figure(1)
% %surface(dif_modos)
% plot(dif_modos(:,1),'r');
% hold on;
% plot(dif_modos(:,2),'b');
% hold on;
% plot(dif_modos(:,3),'g');
% hold on;
% plot(dif_modos(:,4),'y');
% hold on;
% plot(dif_modos(:,5),'k');
% hold off;
% legend('Mode1','Mode2','Mode3','Mode4','Mode5');
% legend ('show');
% title('Diferencia de los modos ponderados con las frecuencias')
%
% Se extrae la parte de interÃ©s (longitud de la viga) de la CWT
%
%
if tam_reduc~=tam_modo;
    for j =1:N_scale
        for i=tam_interpol:tam_interpol*2-1
            coefs_dif_modo1_ext(i-tam_interpol+1,j)=coefs_dif_modo1(j,i);
            coefs_dif_modo2_ext(i-tam_interpol+1,j)=coefs_dif_modo2(j,i);
            coefs_dif_modo3_ext(i-tam_interpol+1,j)=coefs_dif_modo3(j,i);
            coefs_dif_modo4_ext(i-tam_interpol+1,j)=coefs_dif_modo4(j,i);
            coefs_dif_modo5_ext(i-tam_interpol+1,j)=coefs_dif_modo5(j,i);
        end
        for i=tam_interpol:tam_interpol*2-1
        x(i-tam_interpol+1)=x_ext(i);
        end
    end
else
     for j =1:N_scale
        for i=(tam_reduc):(tam_reduc*2-1)
            coefs_dif_modo1_ext(i-tam_reduc+1,j)=coefs_dif_modo1(j,i);
            coefs_dif_modo2_ext(i-tam_reduc+1,j)=coefs_dif_modo2(j,i);
            coefs_dif_modo3_ext(i-tam_reduc+1,j)=coefs_dif_modo3(j,i);
            coefs_dif_modo4_ext(i-tam_reduc+1,j)=coefs_dif_modo4(j,i);
            coefs_dif_modo5_ext(i-tam_reduc+1,j)=coefs_dif_modo5(j,i);
        end
        for i=(tam_reduc):(tam_reduc*2-1)
          x(i-tam_reduc+1)=x_ext(i);
        end
    end
   
end

% 
%
%
% Ploteamos resultados
%
y=1:N_scale;
tam_CWT=length(coefs_dif_modo1_ext(:,1));
L=1.28;
for i=1:length(x)
chi(i)=x(i)/L;
end

 figure(2)
 title('CWT de las diferencias del modo 1 en zona de interï¿½s trï¿½s extensiÃ³n')
 colormap('hot');
 xlabel('Scale');
 ylabel('Position along beam (x/L)');
 ylim([0 1]);
 surface(y,chi,abs(coefs_dif_modo1_ext))
 colorbar;
 print  -r600 -depsc mode1
% 
 figure(3)
 title('CWT de las diferencias del modo 2 en zona de interï¿½s trï¿½s extensiï¿½n')
 colormap('hot');
 xlabel('Scale');
 ylabel('Position along beam (x/L)');
 ylim([0 1]);
 surface(y,chi,abs(coefs_dif_modo2_ext))
 colorbar;
 print  -r600 -depsc mode2
% 
 figure(4)
 title('CWT de las diferencias del modo 3 en zona de interï¿½s trï¿½s extensiï¿½n')
 colormap('hot');
 xlabel('Scale');
 ylabel('Position along beam (x/L)');
 ylim([0 1]);
 surface(y,chi,abs(coefs_dif_modo3_ext))
 colorbar;
 print  -r600 -depsc mode3
% 
 figure(5)
 title('CWT de las diferencias del modo 4 en zona de interï¿½s trï¿½s extensiï¿½n')
 colormap('hot');
 xlabel('Scale');
 ylabel('Position along beam (x/L)');
 ylim([0 1]);
 surface(y,chi,abs(coefs_dif_modo4_ext))
 colorbar;
 print  -r600 -depsc mode4
% 
 figure(6)
 title('CWT de las diferencias del modo 5 en zona de interï¿½s trï¿½s extensiï¿½n')
 colormap('hot');
 xlabel('Scale');
 ylabel('Position along beam (x/L)');
 ylim([0 1]);
 surface(y,chi,abs(coefs_dif_modo5_ext))
 colorbar;
 print  -r600 -depsc mode5

%
% Calculamos la suma de los valores absolutos de las CWT de las diferencias de los modos 1 a 5
% ponderados 
%
for j=1:N_scale
    for i=1:tam_CWT
% Se puede ejecutar tanto el sumatorio de los valores absolutos de los CWT
% de las diferencias de los modos como...
%
        suma_CWT_dif_pond_mod_ext(i,j)=abs(coefs_dif_modo1_ext(i,j))+...
            abs(coefs_dif_modo2_ext(i,j))+abs(coefs_dif_modo3_ext(i,j));%+...
    %      abs(coefs_dif_modo4_ext(i,j))+abs(coefs_dif_modo5_ext(i,j));
    %% ESTA CONFIGURACIï¿½N SIRVE PARA EXALTAR LOS DAï¿½OS CENTRALES
%         suma_CWT_dif_pond_mod_ext(i,j)=abs(coefs_dif_modo1_ext(i,j))+...
%             +abs(coefs_dif_modo3_ext(i,j));%+...
           %abs(coefs_dif_modo4_ext(i,j))+abs(coefs_dif_modo5_ext(i,j));
%
    %% ESTA CONFIGURACIï¿½N SIRVE PARA EXALTAR LOS DAï¿½OS EN BORDE
%         suma_CWT_dif_pond_mod_ext(i,j)=abs(coefs_dif_modo1_ext(i,j))+...
%             +abs(coefs_dif_modo3_ext(i,j));%+...
           %abs(coefs_dif_modo4_ext(i,j))+abs(coefs_dif_modo5_ext(i,j));
%
%           suma_CWT_dif_pond_mod_ext(i,j)=abs(coefs_dif_modo5_ext(i,j))+...
%           abs(coefs_dif_modo4_ext(i,j));
% el valor absoluto del sumatorio de las diferencias de los modos
% ponderados con su signo.
%               suma_CWT_dif_pond_mod_ext(i,j)=abs(coefs_dif_modo1_ext(i,j)+...
%               coefs_dif_modo2_ext(i,j)+coefs_dif_modo3_ext(i,j)+...
%               coefs_dif_modo4_ext(i,j)+coefs_dif_modo5_ext(i,j));
    end
end

figure(7)
colormap('hot');
xlabel('Scale','FontSize',14);
ylabel('Position along beam (x/L)','FontSize',14);
ylim([0 1]);
surface(y,chi,suma_CWT_dif_pond_mod_ext)
colorbar;
print  -r600 -dpng modesum
title('Suma de la CWT de las diferencias de los modos PONDERADOS')
