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

enter image description here

x in blue, y in green, xprime in red

  • $\begingroup$ In the example, $f$ is not linear! As it is $y^2(t)$. Am I right? $\endgroup$ – learner May 25 '14 at 12:03
  • $\begingroup$ @learner: correct, this is a bit harder example. try it either way :) $\endgroup$ – Alex I May 25 '14 at 20:28
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 30 '14 at 10:44
  • $\begingroup$ @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 :) $\endgroup$ – 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))$:


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


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


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;

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

enter image description here

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


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

enter image description here

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;
  • 1
    $\begingroup$ @Downvoter: comments pointing out shortcomings etc. are always appreciated :) $\endgroup$ – Matt L. May 30 '14 at 18:57
  • $\begingroup$ 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)). $\endgroup$ – Alex I May 31 '14 at 1:16
  • $\begingroup$ Btw the downvote wasn't me, obviously I upvoted since this is a great try. $\endgroup$ – Alex I Jun 1 '14 at 2:16
  • $\begingroup$ @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. $\endgroup$ – Matt L. Jun 2 '14 at 9:08
  • $\begingroup$ @AlexI: Please check the update of my answer. $\endgroup$ – Matt L. Jun 2 '14 at 9:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.