# Approximate a polynomial function using a sum of sine waves

I have a polynomial function which I need to approximate by a sum of sine waves with constant amplitude along a given domain.

From what I hear, this might be a good time to make use of Fourier Transforms. However, my function is a simple low-degree polynomial (typically just degree 2) and I need this to be computed very quickly on an everyday computer. I'm not familiar with Fourier Transforms, so I don't know if they're good for this or not.

So here's my approach as of yet: Let's take the function $f(x) = x - \frac16 x^3 + \frac1{180} x^5$. looking around $-\pi$ to $+\pi$, it appears to have an angular frequency of approximately 1 (This is just a slightly modified expansion of the taylor series for sin(x))

Taking advantage of the fact that sin(x) is approximately equal to x around 0, we can estimate f(x) with $\sin(x) - \frac16 \sin(x)^3 + \frac1{180} \sin(x)^5$. But this is not the best estimation. The closer we get to 0, the more sin(x) looks like x. So if we instead use

$L \sin(x/L) - \frac16 (L \sin(x/L))^3 + \frac1{180} (L \sin(x/L))^5$

and make L a large constant, say 1000, the approximation looks very accurate on the given domain. This simplifies to:

$L \sin(x/L) - \frac14\frac{L^3}6(3 \sin(x/L) - \sin(3x/L)) + \frac1{16}\frac{L^5}{180} (10 \sin(x/L) - 5 \sin(3x/L) + \sin(5x/L))$,

or $(L - \frac{L^3}8 + \frac{L^5}{288}) \sin(x/L) + (\frac{L^3}{24} - \frac{L^5}{576}) \sin(3x/L)+\frac{L^5}{2880} \sin(5x/L)$

Now, the problem: This is an extremely good approximation on the domain, but the frequency of these sign waves carries very little meaning because they're mostly arbitrary based on my value of L. There is no unique mapping from the polynomial to the frequencies encoded in it. For instance, I could claim my polynomial is composed of waves of frequencies [0.1, 0.3, and 0.5], or [0.001, 0.003, and 0.005], based on what I set L to. And with sufficiently large L, all the frequencies are effectively 0.

However, most signal processing algorithms all approximately agree with eachother on what frequencies reside in a signal, unlike mine. And it's clear that f(x) is approximately sin(x) on the domain. So it seems like there should be a way to determine that the frequencies in this signal are all around 1.

What are your recommendations for estimating a polynomial function with a sum of sine waves?

As a final note, I did try to simplify my application for the sake of making a concise question. However, if more information helps, my actual function is $g(x) = P(x) \sin(Bx+c)$ where P(x) is a polynomial function and B and c are constants. My goal is to re-express g as a sum of sine waves. If I can express P as a sum of sine waves (which is is what the body of my post here is about), then I can multiply each sine wave by sin(Bx+c) and make use of the fact that $\sin(a) \sin(b) = \frac12 (\cos(a-b)-\cos(a+b))$ to write g as a sum of sine waves. Also, the resultant sine waves do not need to have constant frequencies.

-
You mention in a comment that you don't think Fourier series are what you're looking for, but I do still recommend you test some examples. For your given polynomial $x + x^3/6 + x^5/180$ the first three terms of its Fourier series over the interval $(-\pi,\pi)$ are $$\frac{1}{90} \left(480-50 \pi ^2+\pi ^4\right)\sin(x) + \frac{1}{360} \left(-465+70 \pi ^2-2 \pi ^4\right)\sin(2x) + \frac{1}{7290} \left(5440-870 \pi ^2+27 \pi ^4\right)\sin(3x)$$ which, away from the endpoints at least, gives an okay first approximation. – Antonio Vargas Feb 27 '14 at 15:29
You'll get slightly better results than the one above if you take the endpoints of your interval of interest to be roots of the polynomial. Also, in regards to your idea to use $L\sin(x/L)$, when $L$ is large this is approximately $x$ (to first order), so you're not taking advantage of anything "sinusoidal" by making the substitution. You may as well have replaced $x$ with $L(x/L - (x/L)^3/6)$. It should be noted that there's nothing inherently "sinusoidal" about polynomials--they don't have frequencies to begin with. – Antonio Vargas Feb 27 '14 at 15:37
@Antonio you're right, that is a decent approximation. And the only reason I want a sinusoidal equation is for later post-processing. This is for a form of additive synthesis, where I'll later be doing things like filtering out specific frequencies. I know that amplitude modulation can indirectly affect the frequency of a wave, so my first step was to make all waveforms have constant amplitude. [x] * sin(x) isn't reducible, to my knowledge, but [L*sin(x/L)] * sin(x) is. – Wallacoloo Feb 27 '14 at 17:05
This is a pretty central step in my process, so I'm going to take multiple approaches and see which yields the best results. I'm thinking that I may make two sine waves with identical frequency whose amplitude adds up to the maximum value of my amplitude function (say, f(x)=x) over its interval, but initially place them out of phase so they destructively interfere. Then over time I'll bring them back in phase so they constructively interfere and the amplitude will rise. I'm not yet sure if this is possible without creating a phase change in the overall wave though. – Wallacoloo Feb 27 '14 at 17:11
The chosen cutoff point is all essential. When you have a signal that is essentially $\equiv0$ for $|x|>L$ then all cutoffs at $L'\gg L$ produce similar results; but this is not the case for your polynomial. – Christian Blatter Mar 6 '14 at 21:06

Representing a signal as a sum of sine waves is exactly what the Fourier transform does and it can be computed extremely fast with the Fast Fourier Transform (FFT) algorithm.

Evaluate your function on a set of $N$ evenly spaced points. For instance, let $x_n = -\pi + n \Delta x$ for $n=0,1,\dots,N-1$ and $\Delta x = 2\pi / (N-1)$ similar to your example. Then take the the FFT of $f(\mathbf{x})$. The result of will be the weights of a set of complex exponentials such that their weigthed sum is equal to $f(\mathbf{x})$. That is, if $\mathbf{w} = \operatorname{FFT}(f(\mathbf{x}))$, then

$$f(x_n) = \sum_{m=0}^{N-1} w_m e^{j 2\pi n m/N}$$

If you want to stay away from complex exponentials, you can use the Discrete Cosine Transform instead, which as the name suggests will only use cosines in the weighted sum and will look something like

$$f(x_n) = \sum_{m=0}^{N-1} w_m \cos \left( {\pi\over N}\left(m+{1\over 2}\right) n \right).$$

-
Thanks. I actually just looked up Fourier transforms last night, and they were much easier to understand than I would have expected. It appears that what this transform does that's significantly different from my approach is that it assumes that the signal is cyclic and that the length of the domain is the length of one cycle. This indirectly adds a requirement that no frequency can have an amplitude larger than the amplitude of the overall signal on the domain (correct me if I'm wrong) - something which my approach violates. – Wallacoloo Feb 27 '14 at 0:56
And so I think there truly are an infinite number of ways to model a function along a finite domain using sine waves. It's only when we introduce the above constraints that the solution becomes unique (I think). Yet signal filters in hardware are able to identify frequencies before they see the entire signal. That is, they're able to do it without the same constraints as a Fourier transform. So I'm going to look into how they are able to produce a unique result without those constraints. – Wallacoloo Feb 27 '14 at 1:01

On the domain of $[-1, 1], x = \cos(\arccos(x))$

We can substitude $x \over b$ in place of the original x in order to make a similar equality, $x = b\cos(\arccos{x \over b})$, valid on [-b, b].

This can be generalized to $x^p = b^p\cos(\arccos{x^p \over b^p})$

Therefore, my original function, $f(x) = x - \frac16 x^3 + \frac1{180} x^5$, is equivalent to:

$$f(x) = b\cos(\arccos{x \over b}) - \frac16 b^3\cos( \arccos{x^3 \over b^3}) + \frac1{180} b^5 \cos(\arccos{x^5 \over b^5})$$

along the domain, [-b, b].

It's a blindingly obvious solution to me, and it appears to fulfill all of my (perhaps poorly-stated) requirements: the amplitude of each sinewave is constant and even though the frequency varies, this is acceptable because it's still trivial to calculate the instantaneous frequency of each wave at any time. Furthermore, the amplitude of each wave never exceeds that of the original signal, which appeared to be one of the biggest problems in my original approach. My only complaints are that the frequency is undefined at $\pm b$, and the frequencies can vary based on the domain that you look at, though that is a problem with the Fourier transform approach as well.

-