Descargar Audio del Diapasón (ZIP)

Práctica 5: Sensores Inteligentes

Analogías Dinámicas

Objetivo:
  1. Comprender el comportamiento dinámico de un sistema mecánico vibratorio (diapasón) y su analogía con sistemas de segundo orden.
  2. Realizar la adquisición de señales acústicas utilizando herramientas de hardware estándar (micrófono) y MATLAB.
  3. Aplicar técnicas de análisis espectral (FFT) para identificar la frecuencia fundamental de resonancia.
  4. Obtener un modelo matemático de caja negra (Black-box) del sistema mediante algoritmos de identificación de subespacio (N4SID) y validarlo con datos experimentales.
Motivación:
La presente práctica se motiva en la necesidad de caracterizar experimentalmente un oscilador armónico de alto factor de calidad. El diapasón sirve como plataforma experimental para validar modelos matemáticos de sistemas de segundo orden, comprender los principios físicos detrás de la sensórica inercial moderna (MEMS) y aplicar técnicas de procesamiento digital de señales para la estimación de parámetros dinámicos en el dominio de la frecuencia y el tiempo.
Procedimiento:
A continuación, irás construyendo el sistema paso a paso. Lee la explicación y transcribe el código correspondiente en tu archivo de trabajo

Marco Teórico.

El diapasón

En la figura 1 se muestra un diapasón de uso médico, que es empleado para evaluar la capacidad auditiva de personas adultas. Su construcción se realiza con una barra de acero o aluminio en forma de U, la cual se fija a un vástago de manera que permite vibrar libremente los extremos al ser golpeada. Esta vibración produce un tono cuya frecuencia es fija, y se puede determinar por medio de la ecuación:
(1)
En la ecuación (1), E es el módulo de Young expresado en y ρ es la densidad indicada en , ambas características propias del material. La longitud L y el espesor a, ambos expresados en , son características geométricas del diapasón (figura 2). Por lo tanto, la frecuencia a la cual resuena, está determinada por las características del material empleado y su geometría.
Figura 1. Diapason médico hecho en aluminio
Figura 2. Esquema simplificado de un diapasón

Paso 1. Cálculo de la Frecuencia Teórica.

Primero, calcularemos la frecuencia teórica esperada utilizando las propiedades físicas y geométricas del diapasón. Nota: Ajusta los valores siguientes según tu diapasón real.
L=96.3e-3; % mm
a=5.2e-3; % mm
E = 69e9; % Pa
rho = 2.7e3; % kg/m^3 2.7e-3 g/mm^3
f=(1.875^2/(2*pi*L^2))*sqrt((a^2*E)/(12*rho));
disp(['La frecuencia de vibración calculada para el diapasón es: ', num2str(f), ' [Hz]']);
La frecuencia de vibración calculada para el diapasón es: 457.8524 [Hz]

Paso 2. Adquisición de la Señal Acústica.

Al golpear ligeramente el diapasón contra la palma de la mano, o bien con la rodilla, éste comenzará a vibrar y debería de hacerlo con una frecuencia similar a la calculada previamente. Como el diapasón está rodeado de un medio (aire), esta vibración se debe propagar a través de dicho medio y producirá un sonido; el cual se puede registrar empleando un dispositivo móvil (teléfono celular o tableta) y MATLAB Mobile. Con el archivo de audio registrado mediante este procedimiento, podemos verificar si el diapasón ha producido un sonido cuya frecuencia es similar a la frecuencia calculada previamente. Si quisiera grabar el sonido, puede hacer uso del siguiente código:
% Configuración de grabación
%fs = 44100; % Frecuencia de muestreo [Hz]
%duration = 3; % Duración de la grabación [s]
 
% --- OPCIÓN A: Grabar desde micrófono (Descomentar para usar) ---
% disp('Preparando grabación... Golpea el diapasón ahora.');
% pause(1);
% recObj = audiorecorder(fs, 16, 1);
% disp('Grabando...');
% recordblocking(recObj, duration);
% disp('Fin de grabación.');
% R = getaudiodata(recObj); % 'R' será nuestra señal de respuesta
% ---------------------------------------------------------------------
 
% Pre-procesamiento: Eliminar tendencia y media (Detrend)
%R = detrend(R);
 
% Visualización en el tiempo
%figure;
%plot((0:length(R)-1)/fs, R);
%title('Señal de Audio en el Tiempo');
%xlabel('Tiempo [s]'); ylabel('Amplitud');
%grid on;
nom_arch='D:\Instituto de Ingeniería\Sensores inteligentes\PAPIME 2026\daipason_Al.m4a'
nom_arch = 'D:\Instituto de Ingeniería\Sensores inteligentes\PAPIME 2026\daipason_Al.m4a'
[datos,fs]=audioread(nom_arch);
N=length(datos);
t=0:1/fs:(N-1)*(1/fs);
L=datos(:,1);
R=datos(:,2);
figure('Color', 'w');
subplot(2,1,1)
plot(t, L, 'b')
grid on
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
xlabel('tiempo (s)', 'Color', 'k')
ylabel('amplitud (V)', 'Color', 'k')
title('Lado Izquierdo', 'Color', 'k')
subplot(2,1,2);
plot(t, R, 'r')
grid on
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
xlabel('tiempo (s)', 'Color', 'k')
ylabel('amplitud (V)', 'Color', 'k')
title('Lado Derecho', 'Color', 'k')

Paso 3. Análisis Espectral (FFT)

Como es difícil reconocer la frecuencia del sonido a partir de la gráfica de la señal temporal, se puede obtener una representación alterna mediante la obtención del espectro en frecuencia de la señal empleando el análisis de Fourier [1].
s=R;
S = fft(s)/N;
F = (0:(N/2))*(fs/N);
figure('Color', 'w');
P2 = abs(S/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
plot(F, abs(P1), 'r', 'LineWidth', 2)
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
xlabel('Frecuencia (Hz)', 'Color', 'k')
ylabel('Magnitud', 'Color', 'k')
title('Espectro en frecuencia', 'Color', 'k')
 
axis([0 1000 0 1.1*max(abs(P1))])
grid on
% Encontrar el pico en el espectro
[peakValue, peakIndex] = max(abs(P1));
f_pico= F(peakIndex);
 
% Mostrar el pico encontrado
disp(['La frecuencia del pico en el espectro es: ', num2str(f_pico), ' [Hz]']);
La frecuencia del pico en el espectro es: 458.67 [Hz]
delta_f=f-f_pico;
disp(['La diferencia entre la frecuencia calculada y la frecuencia registrada por el micrófono es: ', num2str(abs(delta_f)), ' [Hz]']);
La diferencia entre la frecuencia calculada y la frecuencia registrada por el micrófono es: 0.81763 [Hz]

Paso 4. Estimación de un Modelo Dinámico

Al tratarse de un sistema cuya respuesta es de naturaleza oscilatoria, suponemos que en primera aproximación, podría representarse por una ecuación difrencial ordinaria de segundo orden:
Utilizaremos el algoritmo N4SID (Numerical algorithms for Subspace State Space System Identification) para obtener un modelo matemático directamente de los datos de audio R [2, 3, 4].
nx=2;
data_id = iddata(R,[],1/fs)
data_id =
Time domain data set with 143356 samples. Sample time: 2.26757e-05 seconds Outputs Unit (if specified) y1 Data Properties
Options = n4sidOptions;
opt.N4weight = 'CVA';
opt.N4Horizon = 'auto';
opt.EstimateCovariance = 'true';
 
% Modelo identificado por N4SID
Modelo = n4sid(data_id,nx,'Form','canonical','Ts',0,Options)
Modelo = Continuous-time identified state-space model: dx/dt = A x(t) + K e(t) y(t) = C x(t) + e(t) A = x1 x2 x1 0 1 x2 -8.449e+06 -114.8 C = x1 x2 y1 1 0 K = y1 x1 4.573e+04 x2 7.003e+08 Parameterization: CANONICAL form with indices: 2. Disturbance component: estimate Number of free coefficients: 4 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using N4SID on time domain data "data_id". Fit to estimation data: 98.31% (prediction focus) FPE: 2.065e-06, MSE: 2.064e-06 Model Properties

Paso 5. Validación del Modelo

Comparamos qué tan bien el modelo matemático puede reproducir (predecir) la señal de audio real que grabamos.
kstep=1;
[R_hat,fit_y,ic]=compare(R,Modelo,kstep);
figure('Color', 'w');
compare(data_id, Modelo, kstep);
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
title('Validación: Datos Reales vs Modelo Estimado', 'Color', 'k');
xlabel('Tiempo', 'Color', 'k');
ylabel('Amplitud', 'Color', 'k');
lgd = findall(gcf, 'Type', 'legend');
set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');
 
grid on;
 
figure('Color', 'w');
subplot(2,1,1)
plot(R, '-r', 'LineWidth', 1.5)
hold on
plot(R_hat, '--b', 'LineWidth', 0.1)
hold off
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
xlabel('Tiempo [s]', 'Color', 'k')
ylabel('Amplitud [V]', 'Color', 'k')
title('Señal de amplitud', 'Color', 'k')
lgd = legend('$\mathbf{y}$', '$\mathbf{\hat{y}}$', 'Interpreter', 'Latex');
set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');
grid on
subplot(2,1,2)
plot(R - R_hat, '-r', 'LineWidth', 1.5)
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k', ...
'GridAlpha', 0.2);
xlabel('Tiempo [s]', 'Color', 'k')
ylabel('Amplitud [V]', 'Color', 'k')
title('Error de estimación', 'Color', 'k')
 
grid on

Actividad

  1. ¿Qué tan cercana fue la frecuencia teórica (Paso 1) a la frecuencia experimental (Paso 3)? Si existe una diferencia grande, ¿a qué factor físico (longitud, material, masa del mango) podrías atribuir el error?
  2. Observa la matriz A del modelo de espacio de estados obtenido en el Paso 4. ¿Tiene la forma canónica de un oscilador armónico?
  3. Utilizando el comando damp(Modelo) en el Paso 5, ¿cuál es el factor de amortiguamientoζ (zeta) estimado? ¿Es un sistema subamortiguado, críticamente amortiguado o sobreamortiguado? ¿Esto coincide con lo que escuchas (el sonido dura mucho o poco)?
  4. Basándote en la ecuación (1), si limaras las puntas del diapasón (reduciendo L), ¿la frecuencia del sonido aumentaría o disminuiría?

Referencias

[1] MathWorks, "fft — Fast Fourier transform," MATLAB Documentation. [Online]. Available: https://www.mathworks.com/help/matlab/ref/fft.html
[2] Apuntes de Sensores Inteligentes — Daan Yael. Maldonado-Uriostigue, D. y Roberto G. Ramírez-Chavarría, Online, 2024. https://sitios.iingen.unam.mx/ApuntesSensoresInteligentes/landing.html
[3] D. Y. Maldonado-Uriostigue, D. Martínez-Gutiérrez and R. G. Ramírez-Chavarría, "Data-driven modeling of DC motor dynamics: a benchmark experiment," en Proceeding of the XX Latin American Control Congress (LACC 2025), Cancún, México, Oct. 15-17, 2025: Springer, en prensa.
[4] MathWorks, "N4SID - Subspace-based state-space model identification," MATLAB System Identification Toolbox Documentation, 2023. [Online]. Available: https://www.mathworks.com/help/ident/ref/n4sid.html.