Análisis espectral de mínimos cuadrados - Least-squares spectral analysis

El análisis espectral de mínimos cuadrados ( LSSA ) es un método para estimar un espectro de frecuencias , basado en un ajuste de mínimos cuadrados de sinusoides a muestras de datos, similar al análisis de Fourier . El análisis de Fourier , el método espectral más utilizado en la ciencia, generalmente aumenta el ruido periódico largo en registros con espacios largos; LSSA mitiga estos problemas. En.

LSSA también se conoce como el método de Vaníček después de Petr Vaníček , y como el método de Lomb (o el periodograma de Lomb ) y el método de Lomb-Scargle (o periodograma de Lomb-Scargle ), basado en las contribuciones de Nicholas R. Lomb y, de forma independiente, Jeffrey D. Scargle. Michael Korenberg y Scott Chen y David Donoho han desarrollado métodos estrechamente relacionados .

Antecedentes históricos

Las estrechas conexiones entre el análisis de Fourier , el periodograma y el ajuste por mínimos cuadrados de sinusoides se conocen desde hace mucho tiempo. Sin embargo, la mayoría de los desarrollos se limitan a conjuntos de datos completos de muestras igualmente espaciadas. En 1963, Freek JM Barning de Mathematisch Centrum , Amsterdam, manejó datos desigualmente espaciados mediante técnicas similares, incluido tanto un análisis de periodograma equivalente a lo que ahora se conoce como el método de Lomb, como el ajuste de mínimos cuadrados de frecuencias seleccionadas de sinusoides determinadas a partir de tales periodogramas. , conectado por un procedimiento que ahora se conoce como seguimiento de emparejamiento con seguimiento de emparejamiento post-backfitting o ortogonal.

Petr Vaníček , un geodesista canadiense de la Universidad de New Brunswick , también propuso el enfoque de búsqueda de emparejamiento, que llamó "análisis espectral sucesivo" y el resultado un "periodograma de mínimos cuadrados", con datos igualmente espaciados y desiguales, en 1969. Generalizó este método para dar cuenta de los componentes sistemáticos más allá de una media simple, como una "tendencia secular lineal (cuadrática, exponencial, ...) predicha de magnitud desconocida", y lo aplicó a una variedad de muestras, en 1971.

El método Vaníček fue luego simplificado en 1976 por Nicholas R. Lomb de la Universidad de Sydney , quien señaló su estrecha conexión con el análisis de periodograma . La definición de periodograma de datos desigualmente espaciados fue posteriormente modificada y analizada por Jeffrey D. Scargle del Centro de Investigación Ames de la NASA , quien demostró que con cambios menores se podría hacer idéntica a la fórmula de mínimos cuadrados de Lomb para ajustar frecuencias sinusoidales individuales.

Scargle afirma que su artículo "no introduce una nueva técnica de detección, sino que estudia la confiabilidad y eficiencia de la detección con la técnica más comúnmente utilizada, el periodograma, en el caso de que los tiempos de observación estén espaciados de manera desigual ", y señala además en referencia al ajuste por mínimos cuadrados de sinusoides en comparación con el análisis de periodograma, que su artículo "establece, aparentemente por primera vez, que (con las modificaciones propuestas) estos dos métodos son exactamente equivalentes".

Press resume el desarrollo de esta manera:

Lomb desarrolló un método completamente diferente de análisis espectral para datos muestreados de manera desigual, uno que mitiga estas dificultades y tiene algunas otras propiedades muy deseables, basado en parte en trabajos anteriores de Barning y Vanicek, y además elaborado por Scargle.

En 1989, Michael J. Korenberg, de la Queen's University en Kingston, Ontario, desarrolló el método de "búsqueda ortogonal rápida" para encontrar más rápidamente una descomposición casi óptima de espectros u otros problemas, similar a la técnica que más tarde se conoció como búsqueda de correspondencia ortogonal. . En 1994, Scott Chen y David Donoho de la Universidad de Stanford desarrollaron el método de "búsqueda de bases" utilizando la minimización de la norma de coeficientes L1 para plantear el problema como un problema de programación lineal , para el cual existen soluciones eficientes.

El método Vaníček

En el método de Vaníček, un conjunto de datos discretos se aproxima mediante una suma ponderada de sinusoides de frecuencias determinadas progresivamente, utilizando una regresión lineal estándar o un ajuste por mínimos cuadrados . Las frecuencias se eligen utilizando un método similar al de Barning, pero yendo más allá en la optimización de la elección de cada nueva frecuencia sucesiva al elegir la frecuencia que minimiza el residuo después del ajuste por mínimos cuadrados (equivalente a la técnica de ajuste ahora conocida como búsqueda de coincidencia con pre backfitting). El número de sinusoides debe ser menor o igual al número de muestras de datos (contando senos y cosenos de la misma frecuencia que sinusoides separados).

Un vector de datos Φ se representa como una suma ponderada de funciones de base sinusoidal, tabulada en una matriz A mediante la evaluación de cada función en los tiempos de muestra, con el vector de peso x :

${\ Displaystyle \ phi \ approx {\ textbf {A}} x}$

donde el vector de peso x se elige para minimizar la suma de errores cuadrados al aproximar Φ . La solución para x es de forma cerrada, usando regresión lineal estándar :

${\ Displaystyle x = ({\ textbf {A}} ^ {\ mathrm {T}} {\ textbf {A}}) ^ {- 1} {\ textbf {A}} ^ {\ mathrm {T}} \ phi.}$

Aquí, la matriz A puede basarse en cualquier conjunto de funciones que sean mutuamente independientes (no necesariamente ortogonales) cuando se evalúen en los tiempos de muestreo; Para el análisis espectral, las funciones utilizadas son típicamente senos y cosenos distribuidos uniformemente en el rango de frecuencia de interés. Si se eligen demasiadas frecuencias en un rango de frecuencia demasiado estrecho, las funciones no serán suficientemente independientes, la matriz estará mal acondicionada y el espectro resultante no será significativo.

Cuando las funciones de base en A son ortogonales (es decir, no correlacionadas, lo que significa que las columnas tienen cero productos punto por pares ), la matriz A ^T A es una matriz diagonal; cuando todas las columnas tienen la misma potencia (suma de cuadrados de elementos), entonces esa matriz es una matriz identidad multiplicada por una constante, por lo que la inversión es trivial. Este último es el caso cuando los tiempos de muestra están igualmente espaciados y los sinusoides se eligen para que sean senos y cosenos igualmente espaciados en pares en el intervalo de frecuencia de 0 a medio ciclo por muestra (espaciados por 1 / N ciclo por muestra, omitiendo el seno fases a 0 y frecuencia máxima donde son idénticamente cero). Este caso particular se conoce como la transformada discreta de Fourier , ligeramente reescrito en términos de datos reales y coeficientes.

${\ Displaystyle x = {\ textbf {A}} ^ {\ mathrm {T}} \ phi}$     (Caso DFT para N muestras y frecuencias igualmente espaciadas, dentro de un factor escalar)

Lomb propuso usar esta simplificación en general, excepto para las correlaciones por pares entre bases de seno y coseno de la misma frecuencia, ya que las correlaciones entre pares de sinusoides son a menudo pequeñas, al menos cuando no están muy cerca. Esta es esencialmente la formulación de periodograma tradicional , pero ahora se ha adoptado para su uso con muestras espaciadas de manera desigual. El vector x es una buena estimación de un espectro subyacente, pero dado que se ignoran las correlaciones, A x ya no es una buena aproximación a la señal, y el método ya no es un método de mínimos cuadrados; sin embargo, se ha seguido mencionando como tal.

El periodograma Lomb-Scargle

En lugar de simplemente tomar productos de puntos de los datos con formas de onda seno y coseno directamente, Scargle modificó la fórmula de periodograma estándar para encontrar primero un retardo de tiempo tal que este par de sinusoides serían mutuamente ortogonales en tiempos de muestreo , y también se ajustaron para potencias potencialmente desiguales. de estas dos funciones base, para obtener una mejor estimación de la potencia a una frecuencia, lo que hizo que su método de periodograma modificado sea exactamente equivalente al método de mínimos cuadrados de Lomb. El tiempo de retraso está definido por la fórmula ${\ Displaystyle \ tau}$${\ Displaystyle t_ {j}}$${\ Displaystyle \ tau}$

${\ Displaystyle \ tan {2 \ omega \ tau} = {\ frac {\ sum _ {j} \ sin 2 \ omega t_ {j}} {\ sum _ {j} \ cos 2 \ omega t_ {j}} }.}$

El periodograma en frecuencia se estima entonces como: ${\ Displaystyle \ omega}$

${\ Displaystyle P_ {x} (\ omega) = {\ frac {1} {2}} \ left ({\ frac {\ left [\ sum _ {j} X_ {j} \ cos \ omega (t_ {j } - \ tau) \ right] ^ {2}} {\ sum _ {j} \ cos ^ {2} \ omega (t_ {j} - \ tau)}} + {\ frac {\ left [\ sum _ {j} X_ {j} \ sin \ omega (t_ {j} - \ tau) \ right] ^ {2}} {\ sum _ {j} \ sin ^ {2} \ omega (t_ {j} - \ tau)}} \ right)}$

que Scargle informa tiene la misma distribución estadística que el periodograma en el caso muestreado uniformemente.

A cualquier frecuencia individual , este método da la misma potencia que un ajuste de mínimos cuadrados a las sinusoides de esa frecuencia, de la forma ${\ Displaystyle \ omega}$

${\ Displaystyle \ phi (t) = A \ sin \ omega t + B \ cos \ omega t.}$

El periodograma generalizado de Lomb-Scargle

El periodograma estándar de Lomb-Scargle es válido para un modelo con media cero. Por lo general, esto se aproxima restando la media de los datos antes de calcular el periodograma. Sin embargo, esta es una suposición inexacta cuando la media del modelo (las sinusoides ajustadas) no es cero. El periodograma Lomb-Scargle generalizado elimina este supuesto y resuelve explícitamente la media. En este caso, la función instalada es

${\ Displaystyle \ phi (t) = A \ sin \ omega t + B \ cos \ omega t + C.}$

El periodograma generalizado de Lomb-Scargle también se ha denominado periodograma de media flotante .

Método de "búsqueda ortogonal rápida" de Korenberg

Michael Korenberg de Queen's University en Kingston, Ontario , desarrolló un método para elegir un conjunto escaso de componentes de un conjunto demasiado completo, como componentes sinusoidales para análisis espectral, llamado búsqueda ortogonal rápida (FOS). Matemáticamente, FOS utiliza una descomposición de Cholesky ligeramente modificada en un proceso de reducción de error cuadrático medio (MSER), implementado como una inversión de matriz dispersa . Al igual que con los otros métodos LSSA, FOS evita la mayor deficiencia del análisis discreto de Fourier y puede lograr identificaciones altamente precisas de periodicidades integradas y sobresale con datos espaciados desigualmente; El método de búsqueda ortogonal rápida también se ha aplicado a otros problemas como la identificación de sistemas no lineales.

Método de "búsqueda de bases" de Chen y Donoho

Chen y Donoho han desarrollado un procedimiento llamado búsqueda de bases para ajustar un conjunto escaso de sinusoides u otras funciones de un conjunto demasiado completo. El método define una solución óptima como aquella que minimiza la norma L1 de los coeficientes, de modo que el problema se puede proyectar como un problema de programación lineal , para el cual se encuentran disponibles métodos de solución eficientes.

Método Chi-cuadrado de Palmer

Palmer ha desarrollado un método para encontrar la función que mejor se ajuste a cualquier número de armónicos elegido, lo que permite más libertad para encontrar funciones armónicas no sinusoidales. Este método es una técnica rápida ( basada en FFT ) para realizar análisis de mínimos cuadrados ponderados en datos espaciados arbitrariamente con errores estándar no uniformes. El código fuente que implementa esta técnica está disponible. Debido a que los datos a menudo no se muestrean en tiempos discretos uniformemente espaciados, este método "cuadrículas" los datos llenando escasamente una matriz de series de tiempo en los tiempos de muestreo. Todos los puntos de la cuadrícula que intervienen reciben un peso estadístico cero, equivalente a tener barras de error infinitas a veces entre muestras.

Aplicaciones

La característica más útil del método LSSA es permitir que los registros incompletos se analicen espectralmente , sin la necesidad de manipular el registro o inventar datos que de otro modo no existirían.

Las magnitudes en el espectro LSSA representan la contribución de una frecuencia o período a la varianza de la serie de tiempo . Generalmente, las magnitudes espectrales definidas de la manera anterior habilitan el régimen de nivel de significancia directo de la salida . Alternativamente, las magnitudes en el espectro de Vaníček también se pueden expresar en dB . Tenga en cuenta que las magnitudes en el espectro de Vaníček siguen una distribución β .

La transformación inversa de la LSSA de Vaníček es posible, como se ve más fácilmente escribiendo la transformación directa como una matriz; la matriz inversa (cuando la matriz no es singular) o pseudo-inversa será entonces una transformación inversa; la inversa coincidirá exactamente con los datos originales si las sinusoides elegidas son mutuamente independientes en los puntos de muestra y su número es igual al número de puntos de datos. No se conoce tal procedimiento inverso para el método de periodograma.

Implementación

El LSSA se puede implementar en menos de una página de código MATLAB . En esencia:

"para calcular el espectro de mínimos cuadrados debemos calcular m valores espectrales ... lo que implica realizar la aproximación de mínimos cuadrados m veces, cada vez para obtener [la potencia espectral] para una frecuencia diferente"

Es decir, para cada frecuencia en un conjunto deseado de frecuencias, las funciones seno y coseno se evalúan en los tiempos correspondientes a las muestras de datos, y se toman y normalizan apropiadamente los productos escalares del vector de datos con los vectores sinusoidales; siguiendo el método conocido como periodograma Lomb / Scargle, se calcula un desplazamiento de tiempo para cada frecuencia para ortogonalizar los componentes seno y coseno antes del producto escalar, como lo describe Craymer; finalmente, se calcula una potencia a partir de esos dos componentes de amplitud . Este mismo proceso implementa una transformada de Fourier discreta cuando los datos están espaciados uniformemente en el tiempo y las frecuencias elegidas corresponden a números enteros de ciclos sobre el registro de datos finito.

Este método trata cada componente sinusoidal de forma independiente o fuera de contexto, aunque no sean ortogonales en los puntos de datos; es el método original de Vaníček. Por el contrario, como explica Craymer, también es posible realizar un ajuste por mínimos cuadrados simultáneo o en contexto completo resolviendo una ecuación matricial, dividiendo la varianza total de los datos entre las frecuencias sinusoidales especificadas. Esta solución de mínimos cuadrados matriciales está disponible de forma nativa en MATLAB como el operador de barra invertida .

Craymer explica que el método simultáneo o en contexto, a diferencia de la versión independiente o fuera de contexto (así como la versión de periodograma debido a Lomb), no puede ajustarse a más componentes (senos y cosenos) que muestras de datos, y además que:

"... también pueden surgir graves repercusiones si las frecuencias seleccionadas dan como resultado que algunos de los componentes de Fourier (funciones trigonométricas) se vuelvan casi linealmente dependientes entre sí, produciendo así una N mal condicionada o casi singular. necesario seleccionar un conjunto diferente de frecuencias a estimar (por ejemplo, frecuencias igualmente espaciadas) o simplemente descuidar las correlaciones en N (es decir, los bloques fuera de la diagonal) y estimar la transformada de mínimos cuadrados inversos por separado para las frecuencias individuales ... "

El método de periodograma de Lomb, por otro lado, puede usar un número arbitrariamente alto o una densidad de componentes de frecuencia, como en un periodograma estándar ; es decir, el dominio de la frecuencia se puede sobremuestrear mediante un factor arbitrario.

En el análisis de Fourier, como la transformada de Fourier o la transformada discreta de Fourier , las sinusoides que se ajustan a los datos son todas mutuamente ortogonales, por lo que no hay distinción entre la simple proyección fuera de contexto basada en el producto escalar en funciones de base versus un ajuste por mínimos cuadrados simultáneos en contexto; es decir, no se requiere inversión de matriz para dividir por mínimos cuadrados la varianza entre sinusoides ortogonales de diferentes frecuencias. Este método se prefiere generalmente por su implementación rápida y eficiente de la transformada de Fourier , cuando se dispone de registros de datos completos con muestras igualmente espaciadas.

Ver también

• Transformada de Fourier discreta no uniforme
• Funciones ortogonales
• SigSpec
• Modelo sinusoidal
• Densidad espectral
• Estimación de densidad espectral , para alternativas competidoras

Referencias

enlaces externos

• Descarga gratuita del paquete LSSA , FORTRAN, el método de análisis espectral de mínimos cuadrados de Vaníček, de la Universidad de New Brunswick .
• La descarga gratuita del paquete LSWAVE , MATLAB, incluye el método de análisis espectral de mínimos cuadrados de Vaníček, de la Encuesta geodésica nacional de EE . UU .
• Descarga gratuita del software LSSA (vía ftp), FORTRAN, método de Vaníček, de Natural Resources Canada .