Desde una perspectiva estadística: Transformada de Fourier frente a regresión con base de Fourier

Question

Estoy tratando de entender si la transformada discreta de Fourier da la misma representación de una curva que una regresión utilizando la base de Fourier. Por ejemplo,

library(fda)
Y=daily$tempav[,1] ## my data
length(Y) ## =365

## create Fourier basis and estimate the coefficients
mybasis=create.fourier.basis(c(0,365),365)  
basisMat=eval.basis(1:365,mybasis)
regcoef=coef(lm(Y~basisMat-1))

## using Fourier transform
fftcoef=fft(Y)

## compare
head(fftcoef)
head(regcoef)

La FFT da un número complejo, mientras que la regresión da un número real.

¿Transmiten la misma información? ¿Existe un mapa uno a uno entre los dos conjuntos de cifras?

(Agradecería que la respuesta se escribiera desde la perspectiva del estadístico en lugar de la del ingeniero. Muchos materiales en línea que puedo encontrar tienen jerga de ingeniería por todas partes, lo que los hace menos apetecibles para mí).

Answer 1

Son iguales. He aquí cómo...

Hacer una regresión

Digamos que ajustas el modelo $$ y_t = \sum_{j=1}^n A_j \cos(2 \pi t [j/N] + \phi_j) $$ donde $t=1,\ldots,N$ y $n = \text{floor}(N/2)$ . Sin embargo, esto no es adecuado para la regresión lineal, por lo que en su lugar se utiliza un poco de trigonometría ( $\cos(a + b) = \cos(a)\cos(b) - \sin(a)\sin(b)$ ) y ajustar el modelo equivalente: $$ y_t = \sum_{j=1}^n \beta_{1,j} \cos(2 \pi t [j/N]) + \beta_{2,j}\sin(2 \pi t [j/N]). $$ Ejecución de la regresión lineal en todas las frecuencias de Fourier $\{j/N : j = 1, \ldots, n \}$ te da un montón ( $2n$ ) de betas: $\{\hat{\beta}_{i,j}\}$ , $i=1,2$ . Para cualquier $j$ , si quisieras calcular el par a mano, podrías usar:

$$ \hat{\beta}_{1,j} = \frac{\sum_{t=1}^N y_t \cos(2 \pi t [j/N])}{\sum_{t=1}^N \cos^2(2 \pi t [j/N])} $$ y $$ \hat{\beta}_{2,j} = \frac{\sum_{t=1}^N y_t \sin(2 \pi t [j/N])}{\sum_{t=1}^N \sin^2(2 \pi t [j/N])}. $$ Se trata de fórmulas de regresión estándar.

Hacer una transformada discreta de Fourier

Cuando se ejecuta una transformada de Fourier, se calcula, para $j=1, \ldots, n$ :

\begin{align*} d(j/N) &= N^{-1/2}\sum_{t=1}^N y_t \exp[- 2 \pi i t [j/N]] \\ &= N^{-1/2}\left(\sum_{t=1}^N y_t\cos(2 \pi t [j/N]) - i \sum_{t=1}^N y_t\sin(2 \pi t [j/N])\right) . \end{align*}

Se trata de un número complejo (observe el $i$ ). Para ver por qué se mantiene esa igualdad, hay que tener en cuenta que $e^{ix} = \cos(x) + i\sin(x)$ , $\cos(-x) = \cos(x)$ y $\sin(-x) = -\sin(x)$ .

Para cada $j$ tomando el cuadrado del conjugado complejo se obtiene " periodograma :"

$$ |d(j/N)|^2 = N^{-1}\left(\sum_{t=1}^N y_t\cos(2 \pi t [j/N])\right)^2 + N^{-1}\left(\sum_{t=1}^N y_t\sin(2 \pi t [j/N])\right)^2. $$ En R, el cálculo de este vector sería I <- abs(fft(Y))^2/length(Y) lo que es un poco raro, porque tienes que escalarlo.

También puede definir el " periodograma a escala " $$ P(j/N) = \left(\frac{2}{N}\sum_{t=1}^N y_t \cos(2 \pi t [j/N]) \right)^2 + \left(\frac{2}{N}\sum_{t=1}^N y_t \sin(2 \pi t [j/N]) \right)^2. $$ Claramente $P(j/N) = \frac{4}{N}|d(j/N)|^2$ . En R sería P <- (4/length(Y))*I[(1:floor(length(Y)/2))] .

La conexión entre ambos

Resulta que la conexión entre la regresión y los dos periodogramas es:

$$ P(j/N) = \hat{\beta}_{1,j}^2 + \hat{\beta}_{2,j}^2. $$ ¿Por qué? Porque la base que elegiste es ortogonal/ortonormal. Usted puede mostrar para cada $j$ que $\sum_{t=1}^N \cos^2(2 \pi t [j/N]) = \sum_{t=1}^N \sin^2(2 \pi t [j/N]) = N/2$ . Introdúcelo en los denominadores de tus fórmulas para los coeficientes de regresión y listo.

Fuente: https://www.amazon.com/Time-Analysis-Its-Applications-Statistics/dp/144197864X

Answer 2

Están estrechamente relacionados. Tu ejemplo no es reproducible porque no incluiste tus datos, así que haré uno nuevo. En primer lugar, vamos a crear una función periódica:

T <- 10
omega <- 2*pi/T
N <- 21
x <- seq(0, T, len = N)
sum_sines_cosines <- function(x, omega){
    sin(omega*x)+2*cos(2*omega*x)+3*sin(4*omega*x)+4*cos(4*omega*x)
}
Yper <- sum_sines_cosines(x, omega)
Yper[N]-Yper[1] # numerically 0

x2 <- seq(0, T, len = 1000)
Yper2 <- sum_sines_cosines(x2, omega)
plot(x2, Yper2, col = "red", type = "l", xlab = "x", ylab = "Y")
points(x, Yper)

Ahora, vamos a crear una base de Fourier para la regresión. Tenga en cuenta que, con $N = 2k+1$ realmente no tiene sentido crear más de $N-2$ funciones de base, es decir $N-3=2(k-1)$ sinusoidales y cosenoidales no constantes, ya que los componentes de frecuencia más alta se aliasan en una malla de este tipo. Por ejemplo, un seno de frecuencia $k\omega$ es indistinguible de un costante (seno): consideremos el caso de $N=3$ es decir, $k=1$ . De todas formas, si quieres volver a comprobarlo, sólo tienes que cambiar N-2 a N en el siguiente fragmento y observe las dos últimas columnas: verá que en realidad son inútiles (y crean problemas para el ajuste, porque la matriz de diseño es ahora singular).

# Fourier Regression with fda
library(fda)
mybasis <- create.fourier.basis(c(0,T),N-2)
basisMat <- eval.basis(x, mybasis)
FDA_regression <- lm(Yper ~ basisMat-1)
FDA_coef <-coef(FDA_regression)
barplot(FDA_coef)

Obsérvese que las frecuencias son exactamente las correctas, pero las amplitudes de los componentes distintos de cero no lo son (1,2,3,4). La razón es que el fda Las funciones de base de Fourier se escalan de una forma extraña: su valor máximo no es 1, como sería para la base de Fourier habitual ${1,\sin{ \omega x},\cos{ \omega x},\dots}$ . No es $\frac{1}{\sqrt \pi}$ tampoco, como lo habría sido para la base ortonormal de Fourier, ${\frac{1}{\sqrt {2\pi}},\frac{\sin{ \omega x}}{\sqrt \pi},\frac{\cos{ \omega x}}{\sqrt \pi},\dots}$ .

# FDA basis has a weird scaling
max(abs(basisMat))
plot(mybasis)

Eso se ve claramente:

1. el valor máximo es inferior a $\frac{1}{\sqrt \pi}$
2. la base de Fourier (truncada a la primera $N-2$ ) contiene una función constante (la línea negra), senos de frecuencia creciente (las curvas que son iguales a 0 en los límites del dominio) y cosenos de frecuencia creciente (las curvas que son iguales a 1 en los límites del dominio), como debe ser

Basta con escalar la base de Fourier dada por fda de forma que se obtenga la base de Fourier habitual, conduce a que los coeficientes de regresión tengan los valores esperados:

basisMat <- basisMat/max(abs(basisMat))
FDA_regression <- lm(Yper ~ basisMat-1)
FDA_coef <-coef(FDA_regression)
barplot(FDA_coef, names.arg = colnames(basisMat), main = "rescaled FDA coefficients")

Probemos fft ahora: tenga en cuenta que como Yper es una secuencia periódica, el último punto no añade realmente ninguna información (la DFT de una secuencia es siempre periódica). Por tanto, podemos descartar el último punto al calcular la FFT. Además, la FFT no es más que un algoritmo numérico rápido para calcular la DFT, y la DFT de una secuencia de números reales o complejos es compleja . Por tanto, lo que realmente queremos son los modulos de los coeficientes de la FFT:

# FFT
fft_coef <- Mod(fft(Yper[1:(N-1)]))*2/(N-1)

Multiplicamos por $\frac{2}{N-1}$ para tener la misma escala que con la base de Fourier ${1,\sin{ \omega x},\cos{ \omega x},\dots}$ . Si no escaláramos, seguiríamos recuperando las frecuencias correctas, pero las amplitudes estarían todas escaladas por el mismo factor con respecto a lo que encontramos antes. Tracemos ahora los coeficientes fft:

fft_coef <- fft_coef[1:((N-1)/2)]
terms <- paste0("exp",seq(0,(N-1)/2-1))
barplot(fft_coef, names.arg = terms, main = "FFT coefficients")

Vale: las frecuencias son correctas, pero ten en cuenta que ahora las funciones base ya no son senos y cosenos (son exponenciales complejas $\exp{ni\omega x}$ donde con $i$ I denota la unidad imaginaria). Obsérvese también que en lugar de un conjunto de frecuencias distintas de cero (1,2,3,4) como antes, obtenemos un conjunto (1,2,5). La razón es que un término $x_n\exp{ni\omega x}$ en esta expansión de coeficiente complejo (por tanto $x_n$ es complejo) corresponde a dos términos reales $a_n sin(n\omega x)+b_n cos(n\omega x)$ en la expansión de la base trigonométrica, debido a la fórmula de Euler $\exp{ix}=\cos{x}+i\sin{x}$ . El módulo del coeficiente complejo es igual a la suma en cuadratura de los dos coeficientes reales, es decir, $|x_n|=\sqrt{a_n^2+b_n^2}$ . De hecho, $5=\sqrt{3^3+4^2}$ .

Answer 3

Quiero decir que la primera respuesta tiene un error de desbordamiento ya que la segunda fórmula para los coeficientes DFT es infinita cuando j = N/2 en este valor el argumento para el seno es π * t haciendo el valor cero para todos los t's. β^2,j=∑Nt=1ytsin(2πt[j/N])∑Nt=1sin2(2πt[j/N]) Hay otras regresiones de datos DFT en artículos publicados por ahí en google pero no los he analizado correctamente ya que me están confundiendo. Quizás el autor de la primera respuesta pueda ayudarnos ya que estoy enredado con esto. Gracias de antemano