A non-linear state space model where the transitions come via a kernel regression.

nlss.R

library(tuneR);
library(MASS);

rate = 10000;
width = 2500;
zdim = 400;

play_vec = function(v) {
  wnew = Wave(v * 32760/max(abs(v)), right = numeric(0), samp.rate = rate, bit = 16);
	play(wnew, '/Users/rhall/Downloads/playRWave');
}

load_track = function(name) {
	w = readWave(name);
	m = mono(w, "left");
	d = downsample(m, rate);
	rm(w); rm(m);
	x = d@left[1:(floor(length(d@left)/width)*width)];
	rm(d);
	return(t(matrix(data=x, nr=width, nc=length(x) / width)));
}

M = load_track("mindwave.wav")[100:1500,]
print(dim(M));
S = svd(M);

Z = S$u[,1:zdim] %*% diag(S$d[1:zdim]); # latent states.
C = t(S$v[,1:zdim]); # emission matrix.

# cross validate kernel smoothers to find optimal bandwidths.
hv = 2^(1:10);
res = matrix(data = 0, nr = zdim, nc = length(hv));
X = Z[1:(dim(Z)[1]-1),]
y = Z[2:(dim(Z)[1]),]
D = as.matrix(dist(X));
dmax = max(D)
errl = list();
for(i in 1:length(hv)) {
	h = hv[i]
	K = apply(D/dmax, c(1,2), function(a){exp(-h*a^2)})
	errl[[i]] = y
	diag(K) = 0;
	for(j in 1:zdim) {
		yj = y[,j];
		yjh = (K %*% yj) / (K %*% rep(1, length(yj)));
		res[j,i] = sum((yj - yjh)^2) / length(yj);
		errl[[i]][,j] = yjh - yj
	}
}

minidx = apply(res, 1, which.min)
hs = hv[minidx]

# compute error covariance matrix.
err = y
for(j in 1:zdim) {
	err[,j] = errl[[minidx[j]]][,j];
}

Emu = colMeans(err)
EZ = err - rep(1,dim(err)[1]) %*% t(Emu)
EC = t(EZ) %*% EZ / dim(err)[1]


ng = 100;
Zn = matrix(data = 0, nr = ng, nc = zdim)
Zn[1,] = Z[1,]

for(i in 2:ng) {
	dist = t(apply(Z[1:(dim(Z)[1]-1),], 1, function(x,y){sqrt(sum((x-y)^2))/dmax}, Zn[i-1,]))
	znew = rep(0, zdim)
	for(j in 1:zdim) {
		k = exp(-hs[j]*dist^2);
		znew[j] = k %*% Z[2:(dim(Z)[1]),j] / sum(k)
	}
	Zn[i,] = znew + mvrnorm(1, mu = Emu, Sigma = EC)
}

play_vec(as.vector(t(Zn %*% C)))