# 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.

• You can use en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process to find a basis for the signals in the library. Projection of the sensor signal into the basis functions will tell you which library signals are present. – MBaz Apr 12 '16 at 16:38
• Are the signals $A$ and $B$ always at the same offset from each other? Do you know WHEN in time they occur? Is the only thing you're after an estimate of $k$ and $m$ in $x_{\tt measured} = k \cdot A + m \cdot B$ ? – Peter K. Apr 12 '16 at 17:36
• @PeterK. Yes, I am trying to find k and m in that expression, but k and m are discrete. Also, there would be many other possibilities, C, D, E, etc. The x-axis in the plots is actually not time, the sensor has an integration time like a CCD and includes all events that occur within that integration time. – Tony Ruth Apr 12 '16 at 17:58
• @MBaz I think I understand what you are say. I found a new B' which is orthogonal to A, then I calculated Input . A/|A|^2 and Input . B'/|B'|^2 and got numbers close to 1 and 2 respectively. I have a couple questions though: 1. Does this mean that the number of points on the x-axis must be greater than or equal to the number of signals in the library? 2. What is the measure of reliability of identification? How can I get an answer like "None of the above"? – Tony Ruth Apr 12 '16 at 18:28
• @TonyRuth PeterK's answer is good; the technique I suggest is maybe more general but maybe also unnecessarily complex in your case. After projecting the sensor signal on the orthonormal base, you get numbers that indicate how much of each is present on the sensor signal. You don't need more points than signals AFAICT. If the signal is noisy, you're projecting the noise on each basis signal so the numbers you get are noisy too; the reliability will depend on the kind and power of the noise. – MBaz Apr 13 '16 at 1:08

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')

• Multi Hypothesis MF is the way to go :-). – Royi Apr 13 '16 at 6:26
• @Drazick :Baby steps! :-) – Peter K. Apr 13 '16 at 11:22