# How to deconvolve dependent part of signal from independent part?

I have a problem of the following form. There are two signals, x(t) and y(t). The model for the system is such that: $$x(t) = x'(t) + f(y(t))$$ where $f(y(t))$ is a variable offset introduced by interaction of $x$ with $y$. The function $f$ is unknown. As an additional constraint, $x'$ is expected to vary slowly and $y$ varies rapidly (and slowly as well). $x$ has both slow and fast components as a result. $f$ can probably be approximated pretty well by a linear function: it is monotonic increasing and fairly smooth.

The question is how to estimate $f$ and $x'$ given $x$ and $y$?

Sample data:

t = linspace(0, 100, 100)
y = sin(pi * t / 5) + t/50 + 1
x = sin(pi * t/100) + 0.05 * y * y


x in blue, y in green, xprime in red

• In the example, $f$ is not linear! As it is $y^2(t)$. Am I right? – learner May 25 '14 at 12:03
• @learner: correct, this is a bit harder example. try it either way :) – Alex I May 25 '14 at 20:28
• do you have access to both $x(t)$ and $y(t)$? and you want to determine what $x'(t)$ is? except for the non-linearity, a regular old LMS adaptive filter might do this. – robert bristow-johnson May 30 '14 at 10:44
• @robertbristow-johnson: Hi, yes that is correct, have both x(t) and y(t). How to use LMS? Pehaps you could post an answer and collect the bounty :) – Alex I May 31 '14 at 1:17

First of all note that there is a certain ambiguity in the problem formulation because the approximation of $x(t)$ has to be divided into the part represented by $x^{\prime}(t)$ and the part represented by $f(y(t))$, and this division is not unique. The approach I chose is to start with a parametrization of $f(y(t))$, and approximate the given $x(t)$ as well as possible by this function. The remaining approximation error is then taken care of by $x^{\prime}(t)$. I use a least squares criterion for approximating $x(t)$ by $f(y(t))$.

I use the following parametrization for $f(y(t))$:

$$f(y(t))=a_1y^2(t)+a_2y(t)+a_3\tag{1}$$

with unknown parameters $a_i$, which need to be determined given the known functions $x(t)$ and $y(t)$. Of course other choices for $f(t)$ are possible, but let's for illustration purposes use the above parametrization.

First you need to define a grid of values $t_i$ and set up an overdetermined system of equations:

$$x(t_i)=x^{\prime}(t_i)+a_1y^2(t_i)+a_2y(t_i)+a_3,\quad i=0,1,\ldots N-1\tag{2}$$

with some $N$ (which must of course be much greater than 3, the number of unknowns). Define now a $N\times 3$ matrix $\mathbf{Y}$ with rows $Y_i=[y^2(t_i), y(t_i),1]$, a $N\times 1$ vector $\mathbf{x}$ with elements $\mathbf{x}_i=x(t_i)$, and an unknown $3\times 1$ vector $\mathbf{a}$ with elements $a_i$. The optimal least squares solution for $\mathbf{a}$ is given by

$$\mathbf{a}=(\mathbf{Y}^T\mathbf{Y})^{-1}\mathbf{Y}^T\mathbf{x}\tag{3}$$

which in Matlab or Octave is simply implemented by the command

a = Y\x;


The values $x^{\prime}(t_i)$ are then simply given by

$$x^{\prime}(t_i)=x(t_i)-[a_1y^2(t)+a_2y(t_i)+a_3]\tag{4}$$

with the parameters $a_i$ as determined from (3). Note that $x^{\prime}(t)$ as defined by (4) is zero-mean. This is no restriction of generality because the parameter $a_3$ defining the constant part of $f(y(t))$ can as well be added to $x^{\prime}(t)$. This is part of the ambiguity inherent to the problem.

Using your example, the corresponding Matlab/Octave code would be something like

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
a = Y\x;
f = a(1)*y.*y+a(2)*y+a(3)*ones(100,1);
xprime = x - f;
plot(t,x,t,xprime,t,f)


From which you get the following plot ('x' blue, 'xprime' green, 'f' red):

As you can see, 'xprime' is zero mean, but his could easily be changed by adding the constant term determined by $a_3$ to 'xprime' (and subtracting it from 'f').

UPDATE:

So far no constraints have been imposed on $x^{\prime}$. If smoothness of $x^{\prime}$ is desired, this can be accomplished by adding more equations to the overdetermined system (2). These equations can be of the form

$$w\left[a_1(y^2(t_i)-y^2(t_{i-1}))+a_2(y(t_i)-y(t_{i-1}))\right]=w[x(t_i)-x(t_{i-1})],\\\quad i=1,\ldots,N-1\tag{5}$$

with some positive weight factor $w$ determining the relative importance of these additional smoothness constraints. If the matrix $\mathbf{Y}$ and the vector $\mathbf{x}$ are augmented by the values given in (5), the overdetermined system can be solved as before by solving equation (3). I computed the solution for $w=10$. Furthermore I added the constant term $a_3$ to $x^{\prime}$, which is an arbitrary choice (you might as well attribute the constant term to $f(y(t))$). The result is shown below (same functions and colors as in the above plot):

The plot was created in Octave by the following code:

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
w = 10;
y2 = y(2:100); y1 = y(1:99);
x2 = x(2:100) - x(1:99);
Y2 = w * [ y2.*y2 - y1.*y1, y2 - y1, zeros(99,1)];
Y = [Y;Y2];
X = [x;w*x2];
a = Y\X
f = a(1)*y.*y+a(2)*y;
xprime = x - f;
plot(t,x,t,xprime,t,f+a(3))

• @Downvoter: comments pointing out shortcomings etc. are always appreciated :) – Matt L. May 30 '14 at 18:57
• Thank you. Unfortunately this doesn't recover the real x' even in a toy problem with no noise, so I think it is trying to minimize the wrong thing. The constraint is not that x' is as small as possible (which is what I think you're doing with least squares) but rather that it be varying as slowly as possible. Thus, perhaps the LSE could be on sum(abs(x(t) - x(t-1))) instead of sum(abs(x(t)). – Alex I May 31 '14 at 1:16
• Btw the downvote wasn't me, obviously I upvoted since this is a great try. – Alex I Jun 1 '14 at 2:16
• @AlexI: I wasn't worried about that. I'll add an update to my answer soon. I still believe that the basic approach is appropriate. If you want smoothness constraints on xp then this can be added to the current method. – Matt L. Jun 2 '14 at 9:08
• @AlexI: Please check the update of my answer. – Matt L. Jun 2 '14 at 9:22