Algorithm for Discrete Signal Decomposition 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.
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 &lt;- 40

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


B &lt;- rep(0, prototype_length)
B[(prototype_length/2+1):(prototype_length)] &lt;- `^`(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 &lt;- sample(0:10, 1)
m &lt;- sample(0:10, 1)

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

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

image(0:10,0:10,err, xlab="k hat", ylab = "m hat")
ix &lt;- 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')