Nothing
## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 7,
fig.height = 5
)
## ---- include=FALSE------------------------------------------------------
if (as.character(Sys.info()[8]) == "lee") {
devtools::load_all("/home/lee/Dropbox/swdft/r/swdft")
} else {
library(swdft)
}
## ---- eval=FALSE---------------------------------------------------------
# library(swdft)
## ------------------------------------------------------------------------
N <- 100
window_size <- 32
x_random <- rnorm(n=N, mean=0, sd=1)
a <- swdft::swdft(x=x_random, n=window_size)
## ------------------------------------------------------------------------
# --- Generate a local cosine signal + Gaussian noise ---
set.seed(999)
## Specify the parameters for the local cosine signal
signal_length <- 96 # The length of the local cosine signal
freq <- 3 / 32 # Frequency, interpreted as "3 complete cycles every 32 data-points"
amplitude <- 1 # The amplitude of the local cosine signal
phase <- 0 # Phase of the local cosine signal
periodic_start <- 20 # When the local cosine signal starts
periodic_length <- 50 # How long the local cosine signal lasts
signal <- swdft::local_signal(N=signal_length, A=amplitude,
Fr=freq, phase=phase, S=periodic_start,
L=periodic_length)
## Generate the Gaussian noise factor, and add to the local cosine signal
sigma <- .3 # Standard deviation of the Gaussian noise factor
noise <- rnorm(n=signal_length, mean=0, sd=sigma)
x <- signal + noise
## Plot what the local cosine signal + noise looks like
plot(x, pch=19, cex=1.4, xlab="", ylab="Signal",
main="Local Cosine Signal plus Gaussian Noise")
lines(signal, lwd=2, col="red")
## ------------------------------------------------------------------------
window_size <- 32
a <- swdft::swdft(x=x, n=window_size)
## ------------------------------------------------------------------------
length(x)
dim(a$a)
## ------------------------------------------------------------------------
a_nopad <- swdft::swdft(x=x, n=window_size, pad=FALSE)
dim(a_nopad$a)
## ------------------------------------------------------------------------
amod <- Mod(a$a)^2
## ---- eval=FALSE---------------------------------------------------------
# swdft_fft <- function(x, n, taper) {
# N <- length(x)
# P <- N - n + 1
# a <- array(data = NA, dim = c(n, P))
#
# for (p in n:N) {
# a[, p - n + 1] <- stats::fft(z = x[(p - n + 1):p] * taper)
# }
#
# return(a)
# }
## ------------------------------------------------------------------------
plot(a)
## ------------------------------------------------------------------------
plot(a, col="tim.colors", title="SWDFT using the 'tim.colors' colorscale")
## ------------------------------------------------------------------------
sampling_frequency <- 100
plot(a, freq_type="hertz", fs=sampling_frequency)
## ------------------------------------------------------------------------
plot(a, freq_type="fraction")
## ------------------------------------------------------------------------
years <- 1900:1995
plot(a, freq_type="fraction", col="tim.colors", custom_xaxis=years, xlab="Years")
## ------------------------------------------------------------------------
a_taper <- swdft::swdft(x=x, n=window_size, taper_type="cosine", p=.2)
plot(a_taper)
## ------------------------------------------------------------------------
a_smooth <- swdft::swdft(x=x, n=window_size, smooth="daniell", m=1, num_convs=0)
plot(a_smooth)
## ------------------------------------------------------------------------
cosreg_fit <- swdft::cosreg(x=x, f=freq)
plot(cosreg_fit)
## ------------------------------------------------------------------------
coefficients(cosreg_fit)
## ------------------------------------------------------------------------
periodogram <- Mod(fft(x))^2
freqs <- (0:(length(x) - 1)) / length(x)
max_freq <- freqs[which.max(periodogram)]
cat("Estimated Frequency: ", max_freq, " True Frequency: ", freq, " \n")
## ------------------------------------------------------------------------
swdft::get_max_freq(x=x) == max_freq
## ------------------------------------------------------------------------
local_cosreg_fit <- swdft::local_cosreg(x=x)
plot(local_cosreg_fit)
coefficients(local_cosreg_fit)
## ------------------------------------------------------------------------
set.seed(666)
## Set the frequency and length of the signal
N <- 128
f0 <- 10/N
## Generate a time-varying amplitude
amplitude <- rep(0, N)
inds11 <- 10:20
inds12 <- 21:50
inds13 <- 51:70
amplitude[inds11] <- seq(0, 1, length=length(inds11))
amplitude[inds12] <- seq(1, 1, length=length(inds12))
amplitude[inds13] <- seq(1, 0, length=length(inds13))
## Generate the cosine signal with time-varying amplitude plus noise
signal <- swdft::cosine(N=N, A=amplitude, Fr=f0, phase=0)
noise <- rnorm(n=N, mean=0, sd=sigma)
x_demod <- signal + noise
## Plot what the signal looks like
plot(x_demod, pch=19, cex=1.4, main="Cosine signal with time-varying amplitude plus Gaussian noise",
xlab="", ylab="")
lines(signal, lwd=2, col="red")
## ------------------------------------------------------------------------
complex_demod_fit <- swdft::complex_demod(x=x_demod, f0=f0)
plot(complex_demod_fit)
lines(signal, col="blue")
legend("topright", col=c("black", "red", "blue"),
c("Signal + Noise", "Complex Demodulation", "True Signal"), lwd=1)
## ------------------------------------------------------------------------
complex_demod_fit_smoother <- swdft::complex_demod(x=x_demod, f0=f0, passfreq=.05)
plot(complex_demod_fit_smoother)
lines(signal, col="blue")
legend("topright", col=c("black", "red", "blue"),
c("Signal + Noise", "Complex Demodulation", "True Signal"), lwd=1)
## ------------------------------------------------------------------------
## Generate the time-varying amplitude for the second periodic component
amplitude2 <- rep(0, N)
inds21 <- 50:70
inds22 <- 71:100
inds23 <- 101:120
amplitude2[inds21] <- seq(0, 1, length=length(inds21))
amplitude2[inds22] <- seq(1, 1, length=length(inds22))
amplitude2[inds23] <- seq(1, 0, length=length(inds23))
## Set the frequency for the second periodic component
f1 <- 30 / N
## Generate the second signal and add it to the first
signal2 <- swdft::cosine(N=N, A=amplitude2, Fr=f1, phase=phase)
x_demod <- x_demod + signal2
## Plot what the two signals plus noise look like
plot(x_demod, pch=19, cex=1, main="Two cosine signals with time-varying amplitudes plus noise",
xlab="")
lines(signal, col="red", lwd=2)
lines(signal2, col="blue", lwd=2)
## ------------------------------------------------------------------------
matching_demod_fit <- swdft::matching_demod(x=x_demod, n=70)
plot(matching_demod_fit)
## ------------------------------------------------------------------------
plot(x_demod, cex=1, pch=19, ylim=c(-2, 3), xlab="", ylab="",
main="Fitted values from each iteration of matching demodulation")
lines(signal, col="red")
lines(signal2, col="blue")
lines(matching_demod_fit$iterations$iter_fits[1, ], col="green", lwd=2, lty=2)
lines(matching_demod_fit$iterations$iter_fits[2, ], col="purple", lwd=2, lty=2)
legend("topright", c("Signal + Noise", "Signal 1", "Signal 2", "Iteration 1", "Iteration 2"),
col=c("black", "red", "blue", "green", "purple"), lwd=2, cex=.8)
## ------------------------------------------------------------------------
plot(matching_demod_fit$coefficients$inst_amp[1, ], type="l", lwd=2, col="red", ylim=c(0, 1.5),
main="Time-varying amplitude for matching demodulation iterations", xlab="", ylab="")
lines(matching_demod_fit$coefficients$inst_amp[2, ], lwd=2, col="blue")
legend("topright", col=c("red", "blue"), lwd=2,
paste0("Iteration ", 1:2, ", Frequency: ", round(matching_demod_fit$coefficients$f0, digits=3)))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.