La Descomposición de Modo Empírico (EMD) es una técnica poderosa para analizar señales no lineales y no estacionarias, descomponiéndolas en un conjunto de Modos Intrínsecos (IMFs) y un residuo. A continuación, se presenta una implementación en MATLAB y se explican sus componentes.

Función Principal de Análisis EMD

function [imfs, residual, t] = analyze_signal_emd(input_signal, sampling_frequency)
   % Analiza una serie temporal utilizando EMD.
   %
   % Parámetros:
   %   input_signal: Vector de la serie temporal a analizar (1xn).
   %   sampling_frequency: Frecuencia de muestreo en Hz.
   %
   % Retorna:
   %   imfs: Matriz donde cada columna es un IMF.
   %   residual: El componente residual de la descomposición.
   %   t: Vector de tiempo.

   input_signal = input_signal(:); % Asegurar que es un vector columna
   signal_length = length(input_signal);
   t = (0:signal_length-1) / sampling_frequency; % Crear eje de tiempo

   % Realizar la descomposición EMD
   % Se utiliza pchip para una interpolación más suave y se habilita la visualización
   [imfs, residual] = emd(input_signal, 'Interpolation', 'pchip', 'Display', 0);
   num_modes = size(imfs, 2);

   % Visualización de los resultados
   figure;
   
   % Señal Original
   ax1 = subplot(num_modes + 2, 1, 1);
   plot(t, input_signal);
   title('Señal Original');
   ylabel('Amplitud');
   grid on;

   % IMF Components
   for i = 1:num_modes
       ax_imf = subplot(num_modes + 2, 1, i + 1);
       plot(t, imfs(:,i));
       title(sprintf('Modo Intrínseco %d (IMF %d)', i, i));
       ylabel('Amplitud');
       grid on;
       
       % Análisis de Frecuencia Instantánea y Espectro de Potencia
       [analytic_signal, ~] = hilbert(imfs(:,i));
       instantaneous_frequency = diff(unwrap(angle(analytic_signal))) / (2 * pi) * sampling_frequency;
       
       % Crear figura separada para análisis de frecuencia e IMF
       figure;
       
       ax_inst_freq = subplot(2, 1, 1);
       plot(t(1:end-1), instantaneous_frequency);
       title(sprintf('Frecuencia Instantánea - IMF %d', i));
       ylabel('Frecuencia (Hz)');
       grid on;
       
       ax_psd = subplot(2, 1, 2);
       pwelch(imfs(:,i), 256, 128, 256, sampling_frequency);
       title(sprintf('Densidad Espectral de Potencia - IMF %d', i));
       xlabel('Frecuencia (Hz)'); % Asegurar etiqueta de frecuencia en el último subplot
       grid on;
   end

   % Componente Residual
   ax_res = subplot(num_modes + 2, 1, num_modes + 2);
   plot(t, residual);
   title('Componente Residual');
   xlabel('Tiempo (s)');
   ylabel('Amplitud');
   grid on;

   % Asegurar que todos los subplots compartan el mismo eje X si es posible
   linkaxes([ax1, ax_imf, ax_res], 'x'); % Enlaza los ejes X de la señal original, IMFs y residual

   % Verificación de Reconstrucción
   reconstructed_signal = sum(imfs, 2) + residual;
   figure;
   plot(t, input_signal, 'Color', [0.5 0.5 0.5], 'DisplayName', 'Señal Original'); % Gris más suave
   hold on;
   plot(t, reconstructed_signal, '--', 'Color', [0.8 0 0], 'DisplayName', 'Señal Reconstruida'); % Rojo punteado
   legend('Location', 'best');
   title('Verificación de Reconstrucción de Señal');
   xlabel('Tiempo (s)');
   ylabel('Amplitud');
   grid on;
   
   % Cálculo del Error Cuadrático Medio Raíz (RMSE)
   reconstruction_error_rmse = sqrt(mean((input_signal - reconstructed_signal).^2));
   fprintf('Error de Reconstrucción (RMSE): %.4f\n', reconstruction_error_rmse);
 end
 

Análisis de Funcionalidad

• Optimización de la Descomposición:

  • Interpolación: Se emplea la interpolación de Hermite cúbica por tramos (pchip) para mitigar los efectos de borde.
  • Criterios de Parada: La detección automática de IMF se basa en criterios establecidos para asegurar la calidad de los modos (diferencia entre puntos etxremos y cruces por cero, y convergencia del promedio de las envolturas a cero).
  • Computación Paralela: Para señales extensas, se puede implementar el procesamiento por segmentos para acelerar el cálculo (requiere modificaciones en la función emd).

• Módulos de Análisis de Señal:

  • Frecuencia Instantánea: Se obtiene mediante la Transformada de Hilbert para caracterizar la dinámica frecuencial de cada IMF.
  • Análisis Espectral: Se utiliza el método de Welch para estimar la Densidad Espectral de Potencia (PSD).
  • Validación de Reconstrucción: Se calcula el RMSE para cuantificar la diferencia entre la señal original y la reconstruida a partir de los IMFs y el residuo.

• Mejoras en Visualización:

  • Diseño de Subgráficos Adaptativo: Muestra dinámicamente todos los componentes IMF.
  • Análisis Interactivo: Permite la inspección de datos en puntos temporales específicos mediante la interacción con los gráficos.
  • Ajuste Dinámico de Parámetros: Posibilidad de modificar los parámetros de descomposición a través de una Interfaz Gráfica de Usuario (GUI).

Ejemplo de Uso

% Generación de una señal de prueba
fs_ejemplo = 1000; % Frecuencia de muestreo (Hz)
t_ejemplo = 0:1/fs_ejemplo:1; % Vector de tiempo de 1 segundo
senal_compuesta = sin(2 * pi * 50 * t_ejemplo) + 0.5 * sin(2 * pi * 120 * t_ejemplo) + 0.3 * randn(size(t_ejemplo));

% Ejecución del análisis EMD
[modos_intr, residuo_final, tiempo_final] = analyze_signal_emd(senal_compuesta, fs_ejemplo);

 

Técnicas de Mejora

• Tratamiento de Efectos de Borde:

  Se puede mejorar la extensión de la señal en los extremos mediante reflexión simétrica para reducir artefactos.

  
  % Ejemplo de extensión simétrica
  longitud_senal = length(senal_compuesta);
  % Extender simétricamente la señal
  senal_extendida = [2*senal_compuesta(1) - senal_compuesta(2:-1:1), senal_compuesta, ...
                   2*senal_compuesta(longitud_senal) - senal_compuesta(longitud_senal-1:-1:max(1, longitud_senal-10))];
  % Re-ejecutar EMD con la señal extendida
  % [imf_ext, residual_ext] = emd(senal_extendida);
  % Recortar los IMFs resultantes para que coincidan con la longitud original
  % imf_procesados = imf_ext(longitud_senal+1:end-longitud_senal, :); 
  
  

• Supresión de Mezcla de Modos (Mode Mixing):

  La Descomposición Empírica de Modo Adaptativo (EEMD) introduce ruido blanco de baja amplitud para promediar los IMFs y mitigar este problema.

  
  % Esquema de mejora EEMD
  num_ensayos = 100; % Número de realizaciones del ruido
  desviacion_estandar_ruido = 0.2; % Amplitud del ruido añadido
  imfs_eemd_acumulados = zeros(length(senal_compuesta), num_ensayos);
  
  for k = 1:num_ensayos
     senal_con_ruido = senal_compuesta + desviacion_estandar_ruido * randn(size(senal_compuesta));
     [imf_temp, ~] = emd(senal_con_ruido);
     % Se asume que el IMF de interés es el último o se selecciona según criterio
     imfs_eemd_acumulados(:, k) = imf_temp(:, end); 
  end
  
  % Promediar los resultados de los ensayos
  imf_promedio_eemd = mean(imfs_eemd_acumulados, 2); 
  % Este imf_promedio_eemd sería un candidato para el IMF final
  
  

• Filtrado Automatizado de Componentes:

  Seleccionar los IMFs más relevantes basándose en criterios como el coeficiente de correlación con la señal original o la energía contenida.

  
  % Método de filtrado por coeficiente de correlación y energía
  umbral_correlacion = 0.2;
  umbral_energia = 0.05;
  
  % Asumiendo que 'modos_intr' contiene los IMFs calculados previamente
  matriz_correlacion = corrcoef(senal_compuesta, modos_intr);
  matriz_correlacion = matriz_correlacion(2:end, 1); % Coeficientes de correlación con la señal original
  
  % Calcular la energía relativa de cada IMF
  energia_total_senal = sum(senal_compuesta.^2);
  energia_imfs = sum(modos_intr.^2, 1);
  energia_relativa_imf = energia_imfs / energia_total_senal;
  
  % Seleccionar IMFs que cumplen ambos umbrales
  indices_seleccionados = find(abs(matriz_correlacion) > umbral_correlacion & energia_relativa_imf > umbral_energia);
  imfs_filtrados = modos_intr(:, indices_seleccionados);
  
  % Reconstrucción con IMFs filtrados (opcional)
  % senal_reconstruida_filtrada = sum(imfs_filtrados, 2) + residuo_final; 
  
  

Interpretación de Resultados

1. Características de los Modos Intrínsecos (IMFs):

   • IMFs de Alta Frecuencia (primeros componentes): Tienden a representar eventos transitorios, ruido o comopnentes de alta frecuencia en la señal.
   • IMFs de Frecuencia Media: Suelen capturar los componentes oscilatorios principales y las periodicidades subyacentes.
   • IMFs de Baja Frecuencia (últimos componentes y residuo): Reflejan la tendencia general o las variaciones a largo plazo de la señal.

2. Escenarios de Aplicación Comunes:

   • Análisis de Vibraciones Mecánicas: Identificación y separación de frecuencias características asociadas a fallos en maquinaria.
   • Señales Biomédicas: Extracción de patrones relevantes, como el complejo QRS en electrocardiogramas (ECG).
   • Series Temporalees Financieras: Detección de patrones de volatilidad y tendencias en mercados.