I'm a programmer, not a mathematician, but I've got a real-world problem I'm trying to solve that's out of my league, and my Google skills so far have failed me.
I have an analog waveform that's been sampled, part of which contains a sine wave, as shown below:
The sample rate of the data and the frequency of the sine wave are both constant, but unrelated. I need to take this sampled waveform, and "upsample" it, IE, generate more samples between the existing sample points. To do this, I've been using Cubic Spline Interpolation. I've taken an existing C++ library from the following location: http://kluge.in-chemnitz.de/opensource/spline/. The formulas are given on the page. I pass in the existing sampled data as control points, with "X" set to the sample number, and "Y" set to the sample value. Given an input value of X, this library generates me the corresponding "Y" value, so passing in an X of 2.5 would generate a value on the generated spline halfway between samples 2 and 3. All of this is working great so far, and it yields a waveform like the one shown above, but I now need to do one extra thing that I'm a bit stuck on.
What I need to be able to do now is "synchronize" with the the sine wave in the sample data. To do that, I want to detect the point at which the sine wave crosses a given threshold, as shown below:
In other words, given a starting value of X, I want to calculate the next X value at which Y equals a given value. Currently I "brute force" an approximate solution, by incrementing by a given step size between two X values until I cross the threshold point. There are two problems with this:
- Accuracy. The solution is only approximate, and I need a highly accurate result. This leads to a very small step size, while still having a margin of error.
- Performance. I have literally billions of these operations to perform, and I want something that can evaluate in a reasonable period of time. Right now I have to trade off accuracy to improve performance.
Instead of using an imperfect brute force search like I currently do, I was wondering if there's a mathematical solution for this problem. I basically want a function that given an X start position and a target Y sample value, will return the next X where Y is equal to the target Y value. I believe this should be possible to achieve, but I lack the math skills to derive the formula myself. Can anyone tell me if such a function is theoretically possible to write, and if so, what a possible implementation would be?
Here are some potentially useful constraints on my input data, which any acceptable solution can assume:
- The X coordinates of control points are linearly increasing values (IE, all evenly spaced and never "doubling back")
- I know the upper and lower bounds for the X and Y coordinates of the control points when generating the spline
Please let me know if you can see a solution to this problem.
Best Answer
Basically, you have a polynomial root-finding problem.
You have a $y$-value, say $y = \bar{y}$, and you want to find $x$ such that $y(x) = \bar{y}$. In other words, you want to find a root of the function $f(x) = y(x) - \bar{y}$. The function $f$ is a cubic polynomial.
A few ideas:
(1) Go find a general-purpose root-finding function. Look in the "Numerical Recipes" book, for example. Brent's method seems to be a popular choice. You can just plug in your function $f$ directly, so less code to write than in the suggestions below.
(2) Go find some polynomial root-finding code. Again, there's some in the "Numerical Recipes" book. This may require you to find the polynomial coefficients of segments of your cubic spline, though. Not difficult, but it's work.
(3) There are formulas for finding the roots of cubic polynomials (see Wikipedia), but coding them correctly is surprisingly difficult. There's good code here and here. This will probably give you the best performance with cubics. Again, it will require you to find the polynomial coefficients of segments of your cubic spline.
(4) Switch to using quadratic splines instead of cubic ones. Then you only have to find roots of quadratic polynomials, not cubics, so the root-finding problem becomes trivial. This will be a pretty big change, but it will pay off if you need high performance.