Bien, finalmente pude implementar el método en Matlab y verificar su funcionamiento. La idea es que si tenemos una señal muy ruidosa, pero si esta es sincrónica a otra, podemos registrar muchas señales (frames) sincronizadas una con otras y luego hacer el promedio punto a punto para eliminar el ruido. Se utiliza muchísimo en análisis de vibraciones, donde la señal del acelerómetro es muy ruidosa, pero el desbalance (que es lo que queremos buscar) está en sincronía con el ángulo de giro del eje. Entonces si podemos tomar 500 vueltas del eje y acomodarlas para que coincidan en tiempo, al promediar punto a punto eliminamos el ruido y aparece la señal buscada, aunque este escondida dentro del ruido.
Aquí dejo un script de matlab por si a alguien le interesa:
% Script para probar el método "Sinchronous Time Averaging"
% Se deben tomar X cantidad de frames, los cuales deben ser sincrónicos
% para luego poder hacer el promedio y eliminar el ruido
clear all;
%Frecuencia de muestreo de las señales
Fs=22050;
%Periodo de muestreo
Ts=1/Fs;
%Muestras por frames. Debe ser potencia de 2 para la FFT
N=1024; % Nos tomará N/Fs seg cada frame -> 1024/22050 = 46mseg
%Frames. A mayor cantidad de frames mayor la eliminación de ruido por
%promediado
Nfrm=200;
%Frecuencias de F1, F2 y de la DPOAE.
%Fueron elegidas como f2/f1=1.2
Fsin1=1291;
Fsin2=1550;
Fsind=2*Fsin1-Fsin2; % 1034 Hz
Af1 = 1000; %66 db SPL
Af2 = 500; %60 db SPL
Afd = 30; %X db
% Vector de tiempo discreto de 1 frame
t=(0:N-1)*Ts;
% Las tres señales que en teoría capta el micrófono.
f1=Af1*sin(2*pi*Fsin1*t);
f2=Af2*sin(2*pi*Fsin2*t);
fd=Afd*sin(2*pi*Fsind*t);
% Sumo las tres señales
signF12d = (fd+f1+f2);
% Creo los Nfrm señales originales agregándoles ruido
matr=zeros(Nfrm, N);
for x=1:Nfrm
matr(x, :) = signF12d + 1000*randn(size(t));
end
%El STA es tomar de todos los frame el 1er elemento y sacar el promedio,
%luego el segundo elemento, luego el terecro, etc. Obteniendo una señal
%nueva de N puntos obtenidos como el promedio de los Nfrm
filtrSign = zeros(1,N);
for y=1:N
tmp=0;
for x=1:Nfrm
tmp = tmp + matr(x, y);
end
filtrSign(y)=tmp/Nfrm;
end
% Calculo el módulo de la FFT de la señal de 1 frame con ruido, normalizada
% al número de muestras
FFT_1frm = abs(fft(matr(1, 1:N),N)/N);
% Calculo el módulo de la FFT de la ftotal, normalizada al número de
% muestras
FFT_filtr = abs(fft(filtrSign,N)/N);
% creo un vector desde 0 hasta Fs/10 incrementando de a Fs/N
f = Fs/10*linspace(0,1,N/10+1);
%Grafico la FFT de 1 frame cualquiera y de la señal filtrada
figure(1);
subplot(1,2,1);
plot(f,2*(FFT_1frm(1:N/10+1)));
subplot(1,2,2);
plot(f,2*(FFT_filtr(1:N/10+1)));
Aquí los resultados:
La primer señal es un frame con ruido. Se ve f1 y f2 que son los tonos emitidos. Pero no se Fdp (La OAE). Al promediar 200 frames, en la segunda imagen se ve que asoma la OAE en 1034 Hz.
Saludos!
