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 dinmicas de la barra SIN dao','Buscar archivo con las frecuencia y amortiguamiento',...
'Salir');
%Cargar archivo con los vectores de frecuencia y amortiguacin
if opcion==1
uiload
else %Salida del programa
break
end
opcion=menu('Cargar archivos con las propiedades dinmicas de la barra SIN dao','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 dinmicas de la barra CON dao','Buscar archivo con las frecuencia y amortiguamiento',...
'Salir');
%Cargar archivo con los vectores de frecuencia y amortiguacin
if opcion==1
uiload
else %Salida del programa
break
end
opcion=menu('Cargar archivos con las propiedades dinmicas de la barra CON dao','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')
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
