Algorithm for Discrete Signal Decomposition

Question

I am building a sensor and I am trying to understand how to process the signal that it generates. The sensor has a library of reference signals. When 'event a' occurs, it produces signal A, when 'event B' occurs it produces signal B. EX:

It is possible for multiple (but few) events to occur, and their signals would add together linearly. I am looking for an algorithm for signal decomposition which takes advantage of the fact that the decomposition is discrete (It can give integer multiples of the library signals only). I also need some measure of the reliability of the identification.

Ex:

Input: (I made this image by taking A + 2 * B + noise)

Example Output:
A + 2 B  -  95% confidence
A + B    -  3% confidence
A        -  1% confidence
None of the Above - 1% confidence    

I've noticed that I can compare two different possibilities by taking an inner product of the normalized proposed solution and input signals. For instance:

Input Signal . (A + 2 B) = 0.999

but

Input Signal . (A + B) = 0.985

So the first solution is a closer match to the input signal. I thought about using a Newton-Raphson to maximize this inner product, but when the actual sensor is in use it could have 1000s of signals in it's library and that would not be realistic.

Answer 1

I'd just take a KIS (keep it simple) approach as a first step.

Define your (unknown) signal as $$ s(t) = k A(t) + m B(t) + n(t)$$ and then just define the error: $$ e(\tilde{k}, \tilde{m}) = \sum_{\forall t} \left| s(t) - \tilde{k} A(t) - \tilde{m} B(t) \right|^2 $$ and then your estimates $\hat{k}$ and $\hat{m}$ are just chosen as: $$ (\hat{k}, \hat{m}) = \arg \min e(\tilde{k}, \tilde{m}) $$

The R code below does this, and the output seems to do the right thing. Here $n(t)$ is a Gaussian noise with zero mean and standard deviation of $0.1$.

The actual $k$ and $m$ values are at the green circle, the estimates are the blue diagonal bar beneath it (in the rightmost plot).

This appears to work for any values of $k$ and $m$, including zero. You can probably extend it for more signals $C(t)$, $D(t)$, etc. too.

---

R Code Below

# 30067

prototype_length <- 40

A <- rep(0, prototype_length)
A[(prototype_length/4+1):(3*prototype_length/4)] <- sin(2*pi*seq(0,prototype_length/2-1)/prototype_length)


B <- rep(0, prototype_length)
B[(prototype_length/2+1):(prototype_length)] <- `^`(sin(2*pi*seq(0,prototype_length/2-1)/prototype_length),12)*3/8

par(mfrow = c(1,3), pty="s")
plot(A, type="l", col="blue", lwd=5)
lines(B, col="red", lwd=5)
title('Original Signals')

k <- sample(0:10, 1)
m <- sample(0:10, 1)

data <- k*A + m*B + 0.1 * rnorm(prototype_length, 0, 0.1)
plot(data, type="l")
title('Scaled & Added Signals')

err <- array(0,c(11,11))
kvals <- array(0,c(11,11))
mvals <- array(0,c(11,11))
for (khat in 0:10)
{
  for (mhat in 0:10)
  {
    err[khat+1,mhat+1] <- sum('^'(abs(data - khat*A - mhat*B), 2))
    kvals[khat+1,mhat+1] <- khat
    mvals[khat+1,mhat+1] <- mhat
  }
}

image(0:10,0:10,err, xlab="k hat", ylab = "m hat")
ix <- which.min(err)
lines(kvals[ix] + c(-0.2,0.2), mvals[ix]+ c(-0.2,0.2), col="blue", lwd = 10)
points(k, m, col="green", lwd = 10)
title('Error, actual values, and minimum')