Created
August 11, 2017 15:39
-
-
Save mpettis/ea1b17899f4d7f5eba0687643859fa5c to your computer and use it in GitHub Desktop.
Fourier DFT FFT reconstruction from important periodogram contributions
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## This is just a quick exploration, answering the question: | |
## How do I reconstruct a signal by picking the highest-power frequency from a periodogram? | |
## Assumes I keep the constant (DC) term. | |
## ----------------------------------------------------------------------------- | |
## Fourier decomposition exploration | |
## ----------------------------------------------------------------------------- | |
x <- 10 * sin(2 * pi * (1:75) / 25) + 300 + rnorm(1:75, 0, 2.5) | |
x | |
fourier_coeff <- fft(x) | |
fourier_coeff | |
iftt <- fft(fourier_coeff, inverse = TRUE) / length(fourier_coeff) | |
## Get the series | |
x <- dat$X1 | |
## get the periodogram | |
xp <- TSA::periodogram(dat$X1, plot = FALSE) | |
#which.max(xp$spec) | |
## Get the fourier transform of x | |
xf <- fft(x) | |
## Make a filtered fft, picking out only the component with the max power from periodogram | |
## Since if periodogram has N entries, the fft has 2N + 1 entries. If I want to keep the constant | |
## term and the biggest periodogram frequency (call it entry K in periodogram), I need to pick out from the fft: | |
## 1 (constant term) | |
## K + 1 (corresponding to K in periodogram) | |
## 2N + 1 - K (corresponding to K in periodogram, which is the equivalent of the c_{-K} fft coefficient) | |
# xffilt <- ifelse(1:length(xf) %in% which.max(xp$spec), xf, 0) | |
# xffilt <- ifelse(Mod(xf) > 100, xf, 0) | |
xffilt <- ifelse(1:length(xf) %in% c(1, which.max(xp$spec) + 1, length(xf) - which.max(xp$spec)), xf, 0) | |
## reconstruct x | |
xr <- Re(fft(xffilt, inverse=TRUE) / length(xffilt)) | |
## plot it | |
plot(x, type="p") | |
par(new=TRUE) | |
plot(xr, type="l") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment