Andrey Akinshin Add prony-2d-noisy figure  over 8 years ago

Commit id: 54e653dfd230ebc0db3fb470ed2fa07140816be0

deletions | additions      

         

\href{https://github.com/AndreyAkinshin/Prony-analysis-Tech-report/blob/master/r/prony-2d-noisy.R}{Prony-2D Noisy problem}        Binary files /dev/null and b/figures/prony-2d-noisy/prony-2d-noisy.png differ          

height = 700\nwidth = 500         

figures/prony-2d-countour3d/prony-2d-countour3d.gif  In_fact_here_you_can__.tex  subsection_Trigonometric_noisy_real_Prony__.tex  figures/prony-2d-noisy/prony-2d-noisy.png           

# History files  .Rhistory  # Saved Data  .RData  # Example code in package build process  *-Ex.R  # R data files from past sessions  .Rdata  # RStudio files  .Rproj.user/           

source("tools-prony-nd.R")  source("tools-math.R")  source("tools-tables.R")  source("tools.R")  n <- 2  mu <- c(1, 1.5)  delta <- 0.05  x <- sort(exp(-2i * pi * delta * mu))  a <- rep(1, n)  m <- prony.moments(x, a)  noise <- generate.noise.seed(52777, 2 * n, 0.1)  reconstruct <- prony.solve(m + noise)  reconstruct.tr <- prony.tr2d.solve(m + noise, delta)  print.table(x, a, m, reconstruct$x, reconstruct$a, m + noise)  print.table(x, a, m, reconstruct.tr$x, reconstruct.tr$a, m + noise)  # Plotting 2D  t <- seq(-5, 5, by = 0.01)  f1 <- gauss(a[1], mu[1], 1)(t)  f2 <- gauss(a[2], mu[2], 1)(t)  f <- f1 + f2  df <- data.frame(t, f1, f2, f) %>% gather("name", "value", 2:4)  plot <- ggplot(df, aes(t, value, group = name, colour = name)) + geom_line()  plot  ggsave("../figures/prony-2d-noisy/prony-2d-noisy.png", plot)           

gauss <- function(a, x0, s2) { function(x) a * exp(-((x-x0)^2)/(2*s2))}  generate.noise <- function(length, mod) {  noise <- runif(length * 2) * mod  noise <- noise[1:length] + noise[(length + 1):(length * 2)] * 1i  noise  }  generate.noise.seed <- function(seed, length, mod) {  set.seed(seed)  generate.noise(length, mod)  }           

clamp <- function(x, a, b) max(a, min(x, b))  prony.moments <- function(x, a) {  n <- length(a)  sapply(0:(2 * n - 1), function(k) sum(a * x ^ k))  }  prony.solve <- function(m) {  n <- length(m) / 2  u.A <- outer(1:n, 1:n, function(i, j) m[i - j + n] * (-1) ^ (j + 1))  u.B <- m[(n + 1):(2 * n)]  u <- solve(u.A, u.B)  x <- sort(polyroot(c(rev(u), 1) * sapply(1:(length(u)+1), function(i) (-1)^i)))  a.A <- t(sapply(0:(n - 1), function(i) x ^ i))  a.B <- m[1:n]  a <- solve(a.A, a.B)  list(x = x, a = a)  }  prony.tr2d.solve <- function(m, delta) {  m0 <- m[1]; m1 <- m[2]; m2 <- m[3]  M0 <- Mod(m0); M1 <- Mod(m1); M2 <- Mod(m2)    Delta.cos <- clamp((2 * M1 ^ 2 - M0 ^ 2 - M2 ^ 2) / (2 * M0 ^ 2 - 2 * M1 ^ 2), -1, 1)  Delta <- acos(Delta.cos)  D <- max(M0 ^ 2 - 4 * (M1 ^ 2 - M0 ^ 2) ^ 2 / (M2 ^ 2 + 3 * M0 ^ 2 - 4 * M1 ^ 2), 0)  a.r <- (M0 - sqrt(D)) / 2  b.r <- (M0 + sqrt(D)) / 2    p <- a.r * Delta.cos + b.r  q <- a.r * sin(Delta)  theta.cos <- clamp((Re(m1) * p + Im(m1) * q) / (p ^ 2 + q ^ 2), -1, 1)  theta <- -acos(theta.cos)  phi <- theta + Delta   mu.r <- phi / (-2 * pi * delta)  nu.r <- theta / (-2 * pi * delta)  x.r <- exp(-2i * pi * delta * mu.r)  y.r <- exp(-2i * pi * delta * nu.r)  list(x = c(x.r, y.r), a = c(a.r, b.r))  }           

library(xtable)  print.table <- function(x1, a1, m1, x2, a2, m2) {  results <- matrix(nrow = 3, ncol = 4 * n)  results[1,] <- c(x1, a1, m1)  results[2,] <- c(x2, a2, m2)  results[3,] <- results[2,] - results[1,]  results <- prettyNum(results, preserve.width = "common", drop0trailing = T, digits = 3, format = "f")  dim(results) <- c(3, 4 * n)  colnames(results) <- c(paste0("x", 1:n), paste0("a", 1:n), paste0("m", 0:(2 * n - 1)))  rownames(results) <- c("S1", "S2", "dS")  xtable(results)  }