Estou tentando entender se a transformada de Fourier discreta fornece a mesma representação de uma curva que uma regressão usando a base de Fourier. Por exemplo,

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)

A FFT fornece um número complexo, enquanto a regressão fornece um número real.

Eles transmitem a mesma informação? Existe um mapa individual entre os dois conjuntos de números?

(Eu apreciaria se a resposta fosse escrita da perspectiva do estatístico em vez da perspectiva do engenheiro. Muitos materiais on-line que posso encontrar têm jargões de engenharia em todo o lugar, o que os torna menos palatáveis ​​para mim.)

Eles são iguais. Aqui está como ...

Fazendo uma regressão

Diga que encaixa no modelo , onde t = 1 , ... , N e n = andar ( N / 2 ) . Porém, isso não é adequado para regressão linear; portanto, você usa alguma trigonometria ( cos ( a + b ) =

$$y_t = \sum_{j=1}^n A_j \cos(2 \pi t [j/N] + \phi_j)$$ $t=1,\ldots,N$$n = \text{floor}(N/2)$ ) e ajuste o modelo equivalente: y t = n ∑ j = 1 β 1 , j cos ( 2 π t [ j / N ] ) + β 2 , j sin ( 2 π t [ j / N ] ) $\cos(a + b) = \cos(a)\cos(b) - \sin(a)\sin(b)$ $$y_t = \sum_{j=1}^n \beta_{1,j} \cos(2 \pi t [j/N]) + \beta_{2,j}\sin(2 \pi t [j/N]).$$ Executando regressão linear em todas as frequências de Fourier dá-lhe um grupo ( 2 N ) de betas: { β i , j } , i = 1 , 2 . Para qualquer j , se você quiser calcular o par manualmente, poderá usar:$\{j/N : j = 1, \ldots, n \}$$2n$$\{\hat{\beta}_{i,j}\}$$i=1,2$$j$

e β 2,j=Σ N t = 1 ytsin(2πt[j/N]

$$\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])}$$ Estas são fórmulas de regressão padrão. $$\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])}.$$

Fazendo uma transformação discreta de Fourier

Quando você executa uma transformação de Fourier, 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*}$$

$i$$e^{ix} = \cos(x) + i\sin(x)$$\cos(-x) = \cos(x)$$\sin(-x) = -\sin(x)$

$j$

$$|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.$$ Em R, calcular esse vetor seria I <- abs(fft(Y))^2/length(Y), o que é meio estranho, porque você precisa escalá-lo.

$$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.$$ $P(j/N) = \frac{4}{N}|d(j/N)|^2$. Em R isso seria P <- (4/length(Y))*I[(1:floor(length(Y)/2))].

A conexão entre os dois

Acontece que a conexão entre a regressão e os dois periodogramas é:

$$P(j/N) = \hat{\beta}_{1,j}^2 + \hat{\beta}_{2,j}^2.$$ Por quê? Porque a base que você escolheu é ortogonal / ortonormal. Você pode mostrar para cada$j$ naquela $\sum_{t=1}^N \cos^2(2 \pi t [j/N]) = \sum_{t=1}^N \sin^2(2 \pi t [j/N]) = N/2$. Conecte isso aos denominadores de suas fórmulas para os coeficientes de regressão e pronto.

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

Eles estão fortemente relacionados. Seu exemplo não é reproduzível porque você não incluiu seus dados, portanto, criarei um novo. Primeiro de tudo, vamos criar uma função 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)

Agora, vamos criar uma base de Fourier para regressão. Note que, com$N = 2k+1$, não faz muito sentido criar mais do que $N-2$ funções básicas, ou seja, $N-3=2(k-1)$senos e cossenos não constantes, porque os componentes de maior frequência são alias nessa grade. Por exemplo, um seno de frequência$k\omega$ é indistinguível de um costante (seno): considere o caso de $N=3$, ou seja, $k=1$. De qualquer forma, se você quiser verificar novamente, basta alterar N-2para Nno trecho abaixo e observar as duas últimas colunas: você verá que elas são realmente inúteis (e criam problemas para o ajuste, porque a matriz de design agora é 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)

Observe que as frequências são exatamente as corretas, mas as amplitudes de componentes diferentes de zero não são (1,2,3,4). O motivo é que as fdafunções da base de Fourier são dimensionadas de maneira estranha: seu valor máximo não é 1, como seria para a base usual de Fourier${1,\sin{ \omega x},\cos{ \omega x},\dots}$. Não é$\frac{1}{\sqrt \pi}$ ou, como teria sido para a 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)

Você vê claramente que:

1. o valor máximo é menor que $\frac{1}{\sqrt \pi}$
2. a base de Fourier (truncada para a primeira $N-2$ termos) contém uma função constante (a linha preta), senos de frequência crescente (as curvas que são iguais a 0 nos limites do domínio) e cossenos de frequência crescente (as curvas que são iguais a 1 nos limites do domínio), como deveria estar

O simples dimensionamento da base de Fourier fornecida por fda, para que a base usual de Fourier seja obtida, leva a coeficientes de regressão com os 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")

Vamos tentar fftagora: observe que, como Yperé uma sequência periódica, o último ponto não adiciona nenhuma informação (a DFT de uma sequência é sempre periódica). Assim, podemos descartar o último ponto ao calcular a FFT. Além disso, a FFT é apenas um algoritmo numérico rápido para calcular a DFT, e a DFT de uma sequência de números reais ou complexos é complexa . Assim, queremos realmente os módulos dos coeficientes da FFT:

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

Nós multiplicamos por $\frac{2}{N-1}$ para ter o mesmo dimensionamento da base de Fourier ${1,\sin{ \omega x},\cos{ \omega x},\dots}$. Se não escalássemos, ainda assim recuperaríamos as frequências corretas, mas as amplitudes seriam todas escaladas pelo mesmo fator em relação ao que encontramos anteriormente. Vamos agora traçar os 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")

Ok: as frequências estão corretas, mas observe que agora as funções básicas não são mais senos e cossenos (são exponenciais complexas $\exp{ni\omega x}$onde com $i$Denomino a unidade imaginária). Observe também que, em vez de um conjunto de frequências diferentes de zero (1,2,3,4) como antes, obtivemos um conjunto (1,2,5). A razão é que um termo$x_n\exp{ni\omega x}$ nesta complexa expansão do coeficiente $x_n$ é complexo) corresponde a dois termos reais $a_n sin(n\omega x)+b_n cos(n\omega x)$ na expansão da base trigonométrica, por causa da fórmula de Euler $\exp{ix}=\cos{x}+i\sin{x}$. O módulo do coeficiente complexo é igual à soma em quadratura dos dois coeficientes reais, ou seja,$|x_n|=\sqrt{a_n^2+b_n^2}$. Na verdade,$5=\sqrt{3^3+4^2}$.