D'un point de vue statistique: transformée de Fourier vs régression avec base de Fourier

13

J'essaie de comprendre si la transformée de Fourier discrète donne la même représentation d'une courbe qu'une régression utilisant la base de Fourier. Par exemple,

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 donne un nombre complexe, tandis que la régression donne un nombre réel.

Transmettent-ils les mêmes informations? Y a-t-il une carte un à un entre les deux ensembles de nombres?

(J'apprécierais que la réponse soit écrite du point de vue du statisticien au lieu du point de vue de l'ingénieur. De nombreux documents en ligne que je peux trouver ont un jargon d'ingénierie partout, ce qui les rend moins agréables au goût.)

qoheleth
la source
Je ne connais pas votre extrait de code, je ne peux donc pas dire si le problème suivant s'y applique. Cependant, la base DFT est généralement définie en termes de fréquences intégrales ("nombre entier"), tandis qu'une "base de fourier" générale pour la régression pourrait utiliser des rapports de fréquence arbitraires (par exemple, y compris les irrationnels, au moins en arithmétique continue). Cela pourrait également être intéressant.
GeoMatt22
Je pense que tout le monde bénéficierait si vous écrivez votre question en termes mathématiques (par opposition aux extraits de code). Quel est le problème de régression que vous résolvez? Quelles sont les fonctions de base de Fourier que vous utilisez? Vous serez surpris de voir comment les réponses à votre question s'amélioreront.
Yair Daon

Réponses:

15

Ce sont les mêmes. Voici comment...

Faire une régression

Supposons que vous correspondiez au modèle t = 1 , , N et n = plancher ( N / 2 ) . Cela ne convient pas à la régression linéaire, cependant, vous utilisez plutôt une trigonométrie ( cos ( a + b ) = cos

yt=j=1nUNEjcos(2πt[j/N]+ϕj)
t=1,,Nn=sol(N/2) ) et ajustez le modèle équivalent: y t = n j = 1 β 1 , j coscos(une+b)=cos(une)cos(b)-péché(une)péché(b) exécution d'une régression linéaire sur toutes les fréquences de Fourier { j / N : j = 1 , , n } vous donne un tas ( 2 n
yt=j=1nβ1,jcos(2πt[j/N])+β2,jpéché(2πt[j/N]).
{j/N:j=1,,n}2n ) de , i = 1 , 2 . Pour tout j , si vous souhaitez calculer la paire à la main, vous pouvez utiliser:{β^je,j}i=1,2j

et β 2,j=

β^1,j=t=1Nytcos(2πt[j/N])t=1Ncos2(2πt[j/N])
Ce sont des formules de régression standard.
β^2,j=t=1Nytsin(2πt[j/N])t=1Nsin2(2πt[j/N]).

Faire une transformation de Fourier discrète

Lorsque vous exécutez une transformée de Fourier, vous calculez, pour :j=1,,n

d(j/N)=N1/2t=1Nytexp[2πit[j/N]]=N1/2(t=1Nytcos(2πt[j/N])it=1Nytsin(2πt[j/N])).

Il s'agit d'un nombre complexe (notez le ). Pour voir pourquoi cette égalité est valable, gardez à l'esprit que e i x = cos ( x ) + i sin ( x ) , cos ( -ieix=cos(x)+isin(x)cos(x)=cos(x)péché(-X)=-péché(X)

j

|d(j/N)|2=N1(t=1Nytcos(2πt[j/N]))2+N1(t=1Nytsin(2πt[j/N]))2.
Dans R, calculer ce vecteur serait I <- abs(fft(Y))^2/length(Y), ce qui est un peu bizarre, car vous devez le mettre à l'échelle.

P(j/N)=(2Nt=1Nytcos(2πt[j/N]))2+(2Nt=1Nytsin(2πt[j/N]))2.
P(j/N)=4N|d(j/N)|2P <- (4/length(Y))*I[(1:floor(length(Y)/2))]

La connexion entre les deux

Il s'avère que le lien entre la régression et les deux périodogrammes est:

P(j/N)=β^1,j2+β^2,j2.
Why? Because the basis you chose is orthogonal/orthonormal. You can show for each j that t=1Ncos2(2πt[j/N])=t=1Nsin2(2πt[j/N])=N/2. Plug that in to the denominators of your formulas for the regression coefficients and voila.

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

Taylor
la source
1
+1 for the answer and the source. It would also be good if you can demonstrate the result with the R objects I posted.
qoheleth
@qoheleth I'll leave that to you. Just be weary of how fft() doesn't scale the way I wrote (I already mentioned this), that I haven't proved anything with intercepts, and that create.fourier.basis() scales the basis functions weirdly.
Taylor
6

They are strongly related. Your example is not reproducible because you didn't include your data, thus I'll make a new one. First of all, let's create a periodic function:

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)

entrez la description de l'image ici

Now, let's create a Fourier basis for regression. Note that, with N=2k+1, it doesn't really make sense to create more than N2 basis functions, i.e., N3=2(k1) non-constant sines and cosines, because higher frequency components are aliased on such a grid. For example, a sine of frequency kω is indistinguishable from a costant (sine): consider the case of N=3, i.e., k=1. Anyway, if you want to double check, just change N-2 to N in the snippet below and look at the last two columns: you'll see that they're actually useless (and they create issues for the fit, because the design matrix is now 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)

entrez la description de l'image ici

Note that the frequencies are exactly the right ones, but the amplitudes of nonzero components are not (1,2,3,4). The reason is that the fda Fourier basis functions are scaled in a weird way: their maximum value is not 1, as it would be for the usual Fourier basis 1,sinωx,cosωx,. It's not 1π either, as it would have been for the orthonormal Fourier basis, 12π,sinωxπ,cosωxπ,.

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

entrez la description de l'image ici

You clearly see that:

  1. the maximum value is less than 1π
  2. the Fourier basis (truncated to the first N2 terms) contains a constant function (the black line), sines of increasing frequency (the curves which are equal to 0 at the domain boundaries) and cosines of increasing frequency (the curves which are equal to 1 at the domain boundaries), as it should be

Simply scaling the Fourier basis given by fda, so that the usual Fourier basis is obtained, leads to regression coefficients having the expected values:

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

entrez la description de l'image ici

Essayons fftmaintenant: notez que puisqu'il Ypers'agit d'une séquence périodique, le dernier point n'ajoute pas vraiment d'informations (la DFT d'une séquence est toujours périodique). Ainsi, nous pouvons ignorer le dernier point lors du calcul de la FFT. De plus, la FFT n'est qu'un algorithme numérique rapide pour calculer la DFT, et la DFT d'une séquence de nombres réels ou complexes est complexe . Ainsi, nous voulons vraiment les modules des coefficients FFT:

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

On multiplie par 2N-1 afin d'avoir la même mise à l'échelle qu'avec la base de Fourier 1,péchéωX,cosωX,. Si nous n'étions pas à l'échelle, nous récupérerions toujours les fréquences correctes, mais les amplitudes seraient toutes mises à l'échelle par le même facteur par rapport à ce que nous avons trouvé auparavant. Voyons maintenant les coefficients 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")

entrez la description de l'image ici

Ok: les fréquences sont correctes, mais notez que maintenant les fonctions de base ne sont plus des sinus et des cosinus (ce sont des exponentielles complexes expnjeωX, où avec jeJe dénote l'unité imaginaire). Notez également qu'au lieu d'un ensemble de fréquences non nulles (1,2,3,4) comme auparavant, nous avons obtenu un ensemble (1,2,5). La raison en est qu'un termeXnexpnjeωX dans cette expansion de coefficient complexe (donc Xn est complexe) correspond à deux termes réels unensjen(nωX)+bncos(nωX) dans l'expansion de la base trigonométrique, en raison de la formule d'Euler expjeX=cosX+jepéchéX. Le module du coefficient complexe est égal à la somme en quadrature des deux coefficients réels, c'est-à-dire|Xn|=unen2+bn2. À vrai dire,5=33+42.

DeltaIV
la source
1
merci DeltaIV, les données dailysont fournies avec le fdapackage.
qoheleth
@qoheleth je ne savais pas. Ce soir, je modifierai ma réponse en utilisant votre jeu de données, et je clarifierai quelques points.
DeltaIV