Aligning time series

gistfile1.r

align2 <- function(x, y)
{
	# FFT
	pad <- rep(0, length(x))
	xp <- c(pad, x, pad)
	yp <- c(pad, y, pad)
	fx <- fft(xp)
	fy <- fft(yp)
	
	# Correlation
	n <- length(fx)
	u <- fy * Conj(fx)
	if(n %% 2 == 2)
		u[n / 2 + 1] <- 0;
  	conv <- abs(fft(u, T));
	i <- which.max(conv) - 1
	
	# Phase shift
	hn <- floor(n/2)
	freqs <- ((1:n + hn - 1) %% n - hn) / n;
	rot <- exp(2i * i * pi * freqs)
	yp <- Re(fft(fy * rot, T) / n)
	yp
}

align <- function(x)
{
	d <- as.matrix(dist(t(x)))
	ref <- x[,i <- which.min(rowSums(d))]
	apply(x, 2, function(y)
	{
		y <- align2(ref, y)
		if(sum(log(cosh(y-ref))) > sum(log(cosh(-y-ref))))
			y <- y * -1
		y
	})
}

pltaln <- function(x, ...)
{
	d <- zapsmall(align(x, ...))
	d[d == 0] <- NA
	az <- apply(is.na(d), 1, all)
	d <- d[!az,]
	rownames(d) <- NULL
	p <- ggplot(melt(d), aes(X1, value, colour = X2)) + geom_line()
	print(p)
}