Calcular e interpretar la frecuencia instantánea

Soy nuevo en el principio de calcular la frecuencia instantánea y se me ocurrieron muchas preguntas al respecto. Los encontrará en una lista de puntos al final de este texto. El texto puede ser un poco largo, discúlpeme por eso, pero realmente intenté trabajar en ese problema por mi cuenta.

Así que estoy interesado en la frecuencia instantánea de una señal de valor real . El cálculo se realiza con la ayuda de una señal analítica , donde es la transformación de Hilbert de .$f(t)$$x(t)$$z(t) = x(t) + j y(t)$$y(t)$$x(t)$

Para calcular las frecuencias instantáneas de la señal analítica seguí el artículo:$z(t)$

El cálculo de la frecuencia instantánea y el ancho de banda instantáneo por Arthur E. Barns a partir de 1992. En este artículo presenta múltiples métodos para calcular la frecuencia instantánea. Escribo todas las fórmulas que propuso (y usé) en un momento.

Para "aprender", jugué con unas señales muy simples y dos un poco más complejas, en MATLAB, y quería obtener sus frecuencias instantáneas.

Fs = 1000;                                            % sampling-rate = 1kHz
t = 0:1/Fs:10-1/Fs;                                    % 10s 'Timevector'
chirp_signal = chirp(t,0,1,2);                         % 10s long chirp-signal, signal 1
added_sinusoid = chirp_signal + sin(2*pi*t*10);        % chirp + sin(10Hz), signal 2
modulated_sinusoid = chirp_signal .* sin(2*pi*t*10);   % chirp * sin(10Hz), signal 3

Las gráficas en el dominio del tiempo de esas tres señales tienen el siguiente aspecto:

Las gráficas de todas las frecuencias instantáneas que obtuve después de aplicar todos los métodos del documento son las siguientes:

Frecuencias instantáneas de la señal de chirp puro: frecuencias instantáneas de la señal de chirp con sinusoide agregado: frecuencias instantáneas de la señal de chirp modulada: tenga en cuenta que en las tres imágenes, el eje y de los gráficos 3 y 4 se amplía, por lo que las amplitudes de esos las señales son muy pequeñas!

La primera posibilidad de pasar de la señal analítica a la frecuencia instantánea es: donde es la fase instantánea. Creo que este es el método más común utilizado hoy en día, al menos en la página web de MATLAB se calcula de esa manera. El código tiene el siguiente aspecto:

$$\begin{equation} f_{2}(t) = \frac{1}{2\pi} \frac{d}{dt}\theta(t) \end{equation}$$ $\theta(t)$

function [instantaneous_frequency] = f2(analytic_signal,Fs)
    factor =  Fs/(2*pi);
    instantaneous_frequency = factor * diff(unwrap(angle(analytic_signal)));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

En el artículo, Graneros ahora sugiere (o más bien dichas compilaciones) otras cuatro formas de calcular frecuencias instantáneas a partir de la señal analítica. También menciona la fórmula superior, pero es la opinión de que no es práctica debido a las ambigüedades en la fase. Supongo que no conocía el unwrap()método, o para ser más precisos, las matemáticas detrás de él. (Yo mismo aprendí sobre ese método hoy, al mirar otros códigos fuente en frecuencias instantáneas)

En su artículo, la fórmula tiene la etiqueta Número (2), por lo tanto, le di a f (t) el índice 2. Todos los demás índices corresponden de la misma manera a sus números en el documento.

Debido a las ambigüedades en la fase, sugiere más bien:

$$\begin{equation} f_{3}(t) = \frac{1}{2\pi}\cdot\frac{\overbrace{x(t)y'(t)}^\text{a}-\overbrace{x'(t)y(t)}^{\text{b}}}{\underbrace{x(t)^2}_{\text{c}}+\underbrace{y(t)^2}_{\text{d}}} \end{equation}$$ Introduje los símbolos "a", "b", "c" y "d", para facilitar un poco la programación:

function [instantaneous_frequency] = f3(analytic_signal,Fs,T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    diff_x = diff(x);
    diff_y = diff(y);
    factor = Fs/(2*pi);
    a = x(2:end).*diff_y;
    b = y(2:end).*diff_x;
    c = x(2:end).^2;
    d = y(2:end).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Luego Barner da tres fórmulas más que denomina "aproximaciones de frecuencia instantánea":

$$\begin{equation} f_{9}(t) = \frac{1}{2\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{x(t)x(t+T)}_{\text{c}}+\underbrace{y(t)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function[instantaneous_frequency] = f9(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(2*pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = x(1:end-T).*x(1+T:end);
    d = y(1:end-T).*y(1+T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append 0 to return-vector to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{11}(t) = \frac{1}{4\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t-T)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t-T)}^{\text{b}}}{\underbrace{x(t-T)x(t+T)}_{\text{c}}+\underbrace{y(t-T)y(t+T)}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = f11(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(4*pi*T);
    a = x(1:end-2*T).*y(1+2*T:end);
    b = x(1+2*T:end).*y(1:end-2*T);
    c = x(1:end-2*T).*x(1+2*T:end);
    d = y(1:end-2*T).*y(1+2*T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [zeros(1,T) instantaneous_frequency zeros(1,T)];
end

$$\begin{equation} f_{14}(t)=\frac{2}{\pi T}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{(x(t)+x(t+T))^2}_{\text{c}}+\underbrace{(y(t)+y(t+T))^2}_{\text{d}}}\right] \end{equation}$$

function [instantaneous_frequency] = formula14(analytic_signal, Fs, T);
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = 2*Fs/(pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = (x(1:end-T)+x(1+T:end)).^2;
    d = (y(1:end-T)+y(1+T:end)).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

En las 3 aproximaciones, las fórmulas T se establecieron en Fs (T = Fs = 1000 = 1s), como se sugiere en el documento.

Ahora mis preguntas son:

• Las fórmulas f2 y f3 devuelven el mismo resultado para la señal de chirp puro. Creo que es bueno, ya que calculan lo mismo. ¡Los tres métodos de aproximación no devuelven lo mismo, ni siquiera algo que esté cerca! ¿Por qué es ese el caso? (Espero que no sea solo un error de programación ...)
• Aunque devuelven lo mismo, especialmente al final de la trama, comienzan a 'moverse' mucho . ¿Cuál es la explicación para eso? Primero pensé en algo como el alias, pero mi frecuencia de muestreo es bastante alta, en comparación con la frecuencia de las señales, por lo que creo que puede excluirse.

• Al menos f2 y f3 parecen funcionar adecuadamente en una señal de chirp puro, pero todos los métodos, incluidos f2 y f3, parecen fallar horriblemente, cuando se trata de más de una frecuencia en la señal. En realidad, tener más de una frecuencia en una señal es siempre el caso. Entonces, ¿cómo se puede obtener la frecuencia instantánea correcta (más o menos)?

  • De hecho, ni siquiera sé qué esperar cuando hay más de una frecuencia presente en la señal. El cálculo devuelve un número para un punto dado en el tiempo, entonces, ¿qué debería hacer cuando, como aquí, hay más frecuencias presentes? ¿Devuelve el promedio de todas las frecuencias o algo así?

• Y mi pregunta probablemente más importante es, ¿cómo se maneja eso en un software real y elaborado? Digamos que quiero saber la frecuencia instantánea de la señal modulada a 1.75 s, y elegí el método f2, de lo que puedo ser 'afortunado' y obtener un número cercano a 6 [Hz] que probablemente sea la respuesta correcta, o yo escojo mis resultados con algunas muestras al lado y de repente obtengo un resultado cableado, demasiado alto, porque desafortunadamente elegí un valor en el pico. ¿Cómo se puede manejar esto? ¿Postprocesándolo con un filtro medio o incluso mejor? Creo que incluso eso podría ser realmente difícil, especialmente en regiones donde hay muchos picos uno al lado del otro.

Y una última pregunta, no tan importante, ¿por qué es que la mayoría de los trabajos que encuentro sobre frecuencias instantáneas son del área de la geografía, especialmente en el cálculo de eventos sismológicos como los terremotos. El artículo de Barne también toma eso como un ejemplo. ¿No es interesante la frecuencia instantánea en muchas áreas?

Hasta ahora, estoy muy agradecido por cada respuesta, especialmente cuando alguien me da consejos sobre cómo implementarlo en un proyecto de software real ;)

Saludos cordiales, Patrick

No es realmente una respuesta, pero puede ser útil: personalmente encontré que el concepto de frecuencia instantánea solo es útil para señales de banda suficientemente estrecha.

Considere el ejemplo simple de dos ondas sinusoidales estables, digamos 100Hz y 934Hz. En este caso, ciertamente puede definir y calcular la frecuencia instantánea (de la forma que desee), pero ¿cuál debería ser el resultado? ¿Qué posible percepción o propiedades puede tener la frecuencia instantánea que diga algo significativo sobre la señal? Aplicar el concepto de frecuencia instantánea a señales que tienen múltiples frecuencias al mismo tiempo, simplemente no tiene mucho sentido.

Es por eso que obtienes resultados decentes para los barridos pero curvas extrañas para el seno Sweep +. También es la razón por la que ves los meneos en la parte alta del barrido. El ancho de banda de la señal es demasiado alto para asignarle un solo número de frecuencia y, por lo tanto, los resultados saltan.

Al menos f2 y f3 parecen funcionar adecuadamente en una señal de chirp puro, pero todos los métodos, incluidos f2 y f3, parecen fallar horriblemente, cuando se trata de más de una frecuencia en la señal. En realidad, tener más de una frecuencia en una señal es siempre el caso. Entonces, ¿cómo se puede obtener la frecuencia instantánea correcta (más o menos)?

Como sugiere Hilmar, el método de transformación de Hilbert (o "Señal analítica") no funciona en banda ancha porque hay más de un componente de frecuencia. puede hacer este método solo para un solo componente sinusoidal.

entonces, con el enfoque de la Señal Analítica, lo que quiere hacer es utilizar esta identidad:

$$\arctan u - \arctan v = \arctan \left( \frac{u-v}{1+uv} \right)$$

sies lo suficientemente pequeño, de lo que puede derivar la fórmula " " de Barner .f $|u-v|$$f_9$

pero debe haber solo una sinusoide variable en el tiempo en el cálculo de la transformación de Hilbert para hacerlo correctamente. y es mejor alinear el componente "en fase" con la salida de la transformación de Hilbert (que se retrasa con un filtro FIR causal). de lo contrario obtendrás basura.

Wow, qué gran pregunta. Voy a responder primero la pregunta no tan importante:

Y una última pregunta, no tan importante, ¿por qué es que la mayoría de los trabajos que encuentro sobre frecuencias instantáneas son del área de la geografía, especialmente en el cálculo de eventos sismológicos como los terremotos. El artículo de Barne también toma eso como un ejemplo. ¿No es interesante la frecuencia instantánea en muchas áreas?

La razón es que el sistema sismográfico "vibroseis" se utiliza en la industria petrolera para realizar estudios sísmicos. Los camiones que he vinculado vibran de aproximadamente 5 Hz a aproximadamente 90 Hz y se pueden hacer para emitir señales de chirrido. Hay mucho dinero en la industria petrolera, y procesar los retornos de estas señales puede ser muy, muy lucrativo. Por lo tanto, muchas personas han pasado muchas horas analizando tales señales, incluso mirando técnicas de frecuencia instantánea.

---

En cuanto a sus preguntas más importantes: en general, hacer diferencias aritméticas y calcular arcotangentes en señales de tiempo discreto es algo malo . Esto se debe a que las estimaciones de frecuencia en tiempo discreto deben calcularse utilizando "aritmética circular" (aritmética vectorial AKA).$^{\rm TM}$

Mira este artículo.

Los mejores enfoques tienden a usar "promedios ponderados por fase" como se implementa aquí . O aquí para un enlace directo a la matlab .

Lamento dar una respuesta un año después del hecho, pero me topé con esta publicación mientras buscaba artículos sobre este mismo tema. Sus preguntas reflejan los desacuerdos e interpretaciones generalizados de la "frecuencia instantánea" que han plagado el campo desde su inicio. Numerosas personas le dirán, como algunas de las respuestas aquí, que IF solo es aplicable a señales de "banda estrecha" o "monocomponente". De hecho, eso no es cierto: a veces el IF obtenido por la transformación de Hilbert se comporta perfectamente bien para las señales de banda ancha y / o "multicomponentes". Una cantidad que se ha propuesto que evita muchas de estas dificultades es la "frecuencia instantánea promedio ponderada (WAIF)", que se puede medir utilizando un espectrograma.

Ver Loughlin en J. Acoust. Soc. Am., Enero de 1999. Otros buenos documentos sobre IF y conceptos erróneos comunes son Picinbono (IEEE Trans. Sig. Proc., Marzo de 1997) y Vakman (IEEE Trans. Sig. Proc., Abril de 1996).