Próbuję zrozumieć, czy dyskretna transformata Fouriera daje taką samą reprezentację krzywej jak regresja z wykorzystaniem zasady Fouriera. Na przykład,

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)

FFT daje liczbę zespoloną, podczas gdy regresja daje liczbę rzeczywistą.

Czy przekazują te same informacje? Czy istnieje mapa jeden do jednego między dwoma zestawami liczb?

(Byłbym wdzięczny, gdyby odpowiedź została napisana z perspektywy statystyka zamiast z perspektywy inżyniera. Wiele materiałów online, które mogę znaleźć, ma żargon techniczny w całym miejscu, co czyni mnie mniej smakowitym.)

Oni są tacy sami. Oto jak...

Robić regresję

Załóżmy, że do modelu gdzie i . Nie jest to jednak odpowiednie dla regresji liniowej, więc zamiast tego używasz trygonometrii ( ) i dopasuj model równoważny: Uruchamianie regresji liniowej na wszystkich częstotliwościach Fouriera daje ci kilka ( ) bet: , . Dla każdego

$$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)$$\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]).$$ $\{j/N : j = 1, \ldots, n \}$$2n$$\{\hat{\beta}_{i,j}\}$$i=1,2$$j$, jeśli chcesz ręcznie obliczyć parę, możesz użyć:

$$\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])}$$ i Są to standardowe formuły regresji. $$\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])}.$$

Wykonywanie dyskretnej transformaty Fouriera

Po uruchomieniu transformacji Fouriera obliczasz dla :$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*}$$

Jest to liczba zespolona (zwróć uwagę na ). Aby zobaczyć, dlaczego ta równość się utrzymuje, należy pamiętać, że , i .$i$$e^{ix} = \cos(x) + i\sin(x)$$\cos(-x) = \cos(x)$$\sin(-x) = -\sin(x)$

Dla każdego biorąc kwadrat złożonej koniugatu daje „ periodogram :”$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.$$ W R obliczenie tego wektora byłoby I <- abs(fft(Y))^2/length(Y)dziwne, ponieważ trzeba go skalować.

Możesz także zdefiniować „ skalowany periodogram ” Wyraźnie . W R byłoby to .

$$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$P <- (4/length(Y))*I[(1:floor(length(Y)/2))]

Połączenie między dwoma

Okazuje się, że związek między regresją a dwoma periodogramami jest:

j ∑ N t

$$P(j/N) = \hat{\beta}_{1,j}^2 + \hat{\beta}_{2,j}^2.$$ Dlaczego? Ponieważ wybrana podstawa jest ortogonalna / ortonormalna. Możesz pokazać dla każdego że . Podłącz to do mianowników formuł dla współczynników regresji i voila.$j$$\sum_{t=1}^N \cos^2(2 \pi t [j/N]) = \sum_{t=1}^N \sin^2(2 \pi t [j/N]) = N/2$

Źródło: https://www.amazon.com/Time-Analysis-Its-Applications-Statistics/dp/144197864X

Są silnie powiązane. Twój przykład nie jest powtarzalny, ponieważ nie podałeś swoich danych, dlatego zrobię nowy. Po pierwsze, stwórzmy funkcję okresową:

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)

Teraz utwórzmy podstawę Fouriera dla regresji. Zauważ, że przy tak naprawdę nie ma sensu tworzyć więcej niż funkcji bazowych, tj. niesynchroniczne sinus i cosinus, ponieważ składniki o wyższej częstotliwości są aliasowane na takiej siatce. Na przykład sinus o częstotliwości jest nie do odróżnienia od costant (sinus): rozważ przypadek , tj. . W każdym razie, jeśli chcesz dwukrotnie sprawdzić, po prostu przejdź do fragmentu kodu poniżej i spójrz na dwie ostatnie kolumny: zobaczysz, że są one faktycznie bezużyteczne (i powodują problemy z dopasowaniem, ponieważ macierz projektowa jest teraz pojedyncza ).$N = 2k+1$$N-2$$N-3=2(k-1)$$k\omega$$N=3$$k=1$N-2N

# 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)

Zauważ, że częstotliwości są dokładnie właściwe, ale amplitudy składników niezerowych nie są (1,2,3,4). Powodem jest to, że fdafunkcje podstawowe Fouriera są skalowane w dziwny sposób: ich maksymalna wartość nie wynosi 1, tak jak w przypadku zwykłej podstawy Fouriera . Nie jest to również , jak byłoby w przypadku ortonormalnej podstawy Fouriera, .${1,\sin{ \omega x},\cos{ \omega x},\dots}$$\frac{1}{\sqrt \pi}$${\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)

Wyraźnie widać, że:

1. maksymalna wartość jest mniejsza niż$\frac{1}{\sqrt \pi}$
2. podstawa Fouriera (skrócona do pierwszych składników ) zawiera stałą funkcję (czarna linia), sinusy o rosnącej częstotliwości (krzywe, które są równe 0 na granicach domen) i cosinus o rosnącej częstotliwości (krzywe, które są równa 1 na granicach domeny), tak jak powinno być$N-2$

Proste skalowanie podanej przez Fouriera podstawy Fouriera fda, aby uzyskać zwykłą podstawę Fouriera, prowadzi do współczynników regresji o oczekiwanych wartościach:

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")

Spróbujmy fftteraz: zwróć uwagę, że ponieważ Yperjest to sekwencja okresowa, ostatni punkt tak naprawdę nie dodaje żadnych informacji (DFT sekwencji jest zawsze okresowa). W ten sposób możemy odrzucić ostatni punkt przy obliczaniu FFT. Ponadto FFT to szybki algorytm numeryczny do obliczania DFT, a DFT ciągu liczb rzeczywistych lub zespolonych jest złożony . Tak więc naprawdę chcemy modułów współczynników FFT:

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

Mnożymy przez , aby uzyskać takie samo skalowanie jak w przypadku podstawy Fouriera . Gdybyśmy nie skalowali, nadal odzyskalibyśmy prawidłowe częstotliwości, ale wszystkie amplitudy byłyby skalowane o ten sam współczynnik w stosunku do tego, co znaleźliśmy wcześniej. Teraz wykreślmy współczynniki fft:$\frac{2}{N-1}$${1,\sin{ \omega x},\cos{ \omega x},\dots}$

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: częstotliwości są poprawne, ale zauważ, że teraz podstawowymi funkcjami nie są już sinus i cosinus (są złożonymi wykładnikami , gdzie oznaczam jednostkę urojoną). Zauważ też, że zamiast zestawu niezerowych częstotliwości (1,2,3,4), jak poprzednio, otrzymaliśmy zestaw (1,2,5). Powodem jest to, że termin w tym złożonym rozwinięciu współczynnika (a więc jest złożony) odpowiada dwóm rzeczywistym warunkom w trygonometryczne rozszerzenie podstawy, ze względu na formułę Eulera . Moduł współczynnika zespolonego jest równy sumie kwadratury dwóch rzeczywistych współczynników, tj.$\exp{ni\omega x}$$i$$x_n\exp{ni\omega x}$$x_n$$a_n sin(n\omega x)+b_n cos(n\omega x)$$\exp{ix}=\cos{x}+i\sin{x}$ 5=$|x_n|=\sqrt{a_n^2+b_n^2}$ . W rzeczywistości .$5=\sqrt{3^3+4^2}$