what is the right algorithm for ... eh, I don't know what it's called

Question

Suppose I have a long and irregular digital signal made up of smaller but irregular signals occurring at various times (and overlapping each other). We will call these shorter signals the "pieces" that make up the larger signal. By "irregular" I mean that it is not a specific frequency or pattern.

Given the long signal I need to find the optimal arrangement of pieces that produce (as closely as possible) the larger signal. I know what the pieces look like but I don't know how many of them exist in the full signal (or how many times any one piece exists in the full signal). What software algorithm would you use to do this optimization? What do I search for on the web to get help on solving this problem?

Answer 1

Here's a stab at it.

This is actually the easier of the deconvolution problems. It is easier in that you may be able to have a unique answer. The harder problem is that you also don't know what the pieces look like. That case is called blind deconvolution. It is a harder problem and is usually iterative and statistical (ML or MAP), and the solution may not be right.

Luckily, your case is easier, but still not so easy because you have multiple pieces :p

I think that it may be commonly called mixture deconvolution?

So let f[t] for t=1,...N be your long signal. Let h1[t]...hn[t] for t=0,1,2,...M be your short signals. Obviously here, N>>M.

So your hypothesis is that:

(1) f[t] = h1[t+a1[1]]+h1[t+a1[2]] + ...
       +h2[t+a2[1]]+h2[t+a2[2]] + ...
       +....
       +hn[t+an[1]]+h2[t+an[2]] + ...

Observe that each row of that equation is actually hj * uj where uj is the sum of shifted Kronecker delta. The * here is convolution.

So now what?

Let Hj be the (maybe transposed depending on how you look at it) Toeplitz matrix generated by hj, then the equation above becomes:

 (2) F = H1 U1 + H2 U2 + ... Hn Un

subject to the constraint that uj[k] must be either 0 or 1.  

where F is the vector [f[0],...F[N]] and Uj is the vector [uj[0],...uj[N]].

So you can rewrite this as:

 (3) F = H * U 

where H = [H1 ... Hn] (horizontal concatenation) and U = [U1; ... ;Un] (vertical concatenation).

H is an Nx(nN) matrix. U is an nN vector.

Ok, so the solution space is finite. It is 2^(nN) in size. So you can try all possible combinations to see which one gives you the lowest ||F - H*U||, but that will take too long.

What you can do is solve equation (3) using pseudo-inverse, multi-linear regression (which uses least square, which comes out to pseudo-inverse), or something like this

Is it possible to solve a non-square under/over constrained matrix using Accelerate/LAPACK?

Then move that solution around within the null space of H to get a solution subject to the constraint that uj[k] must be either 0 or 1.

Alternatively, you can use something like Nelder-Mead or Levenberg-Marquardt to find the minimum of:

  ||F - H U|| + lambda g(U)

where g is a regularization function defined as:

   g(U) = ||U - U*||

where U*[j] = 0 if |U[j]|<|U[j]-1|, else 1

Ok, so I have no idea if this will converge. If not, you have to come up with your own regularizer. It's kinda dumb to use a generalized nonlinear optimizer when you have a set of linear equations.

In reality, you're going to have noise and what not, so it actually may not be a bad idea to use something like MAP and apply the small pieces as prior.