# Fitting a Damped Sine Wave

A damped sine wave is described by

$$x_{(k)} = A \cdot e^{\alpha \cdot k} \cdot cos(\omega \cdot k + p) \tag{1}$$

with frequency $\omega$ , phase p , initial amplitude A and damping constant $\alpha$ . The $x_{(k)}$ are the samples of the function at equally spaced points in time.

With $x_{(k)}$ given, one often has to find the unknown parameters of the function. This can be achieved for instance with nonlinear approximation or with DFT – methods.

I present a method to find the unknown parameters which is based on linear algebra only.

The equation (1) is the solution to the difference equation
$$x_{(k+2)} = 2 \cdot \sqrt b \cdot cos(\omega) \cdot x_{(k+1)} - b \cdot x_{(k)} \tag{2}$$

This difference equation is a linear equation system for the unknown constants b and $2 \cdot \sqrt b \cdot cos(\omega)$.

The overdetermined linear system
$$\left( \begin{matrix} x_{(k+2)} \\ x_{(k+3)} \\ . \\ . \end{matrix} \right) = \left( \begin{matrix} x_{(k+1)} & x_{(k )} \\ x_{(k+2)} & x_{(k+1)} \\ . & .\\ . & . \end{matrix} \right) \cdot \left( \begin{matrix} 2 \cdot \sqrt b \cdot cos(\omega) \\ -b \end{matrix} \right) \tag{3}$$is solved using linear regression analysis.

There is a pitfall in this calculation: For high oversampling rates, i.e. small resulting $\omega$, the numerical results with additional noise are poor. Performance is best near half the Nyquist rate, which results in approximately 4 samples per sine wave. The $x_{(k)}$ can be rearranged accordingly without losing samples or accuracy, such as taking every 8th sample for the $x_{(k)},x_{(k+1)}, ...$ The resultant $\omega$ has to be corrected by this integer factor.

The $b$ is related to $\alpha$ by

$$e^\alpha\ = \sqrt b \tag{4}$$

For the derivation of this relation see the appendix.

So the $\omega$, $\alpha$ and b are known.

The amplitude A and the phase $p$ are still unknown and calculated by
'mixing with a quadrature carrier' :

The system of equations

$$\left( \begin{matrix} x_{(0)} \\ x_{(1)} \\ . \\ . \end{matrix} \right) = \left( \begin{matrix} e^{\alpha \cdot 0} \cdot \sin{(\omega \cdot 0)} & e^{\alpha \cdot 0} \cdot \cos{(\omega \cdot 0)} \\ e^{\alpha \cdot 1} \cdot \sin{(\omega \cdot 1)} & e^{\alpha \cdot 1} \cdot \cos{(\omega \cdot 1)} \\ . & .\\ . & . \end{matrix} \right) \cdot \left( \begin{matrix} u \\ v \end{matrix} \right) \tag{5}$$

is solved for u and v. With the u, v interpreted as the complex number $u+iv$ we obtain the amplitude and phase

$$A=|u+iv|$$
$$p=\sphericalangle{(u-iv)}$$

In real setups the samples $x_{(k)}$ suffer from additional noise. It turns out that the results for frequency $\omega$ and phase $p$ are very stable in noise. That is not the case for amplitude A nor for the damping factor $\alpha$. These parameters are adjusted with a gradient iteration. Amplitude and damping are slightly varied until the best fit to the sampled data is achieved.

I published a corresponding Matlab script in the 'Mathworks file exchange' system. You can find it by searching for 'matlab fit damped sine wave'. Download is free and you do not need a registration. I hope that this algorithm is useful and download numbers will rise.

Cheers
Detlef

Appendix

Replacing $x_{(k)}$ in (2) with its definition in (1) yields
\begin{align} & A e^{\alpha \cdot (k+2)} \cos {(\omega (k+2) + p)} = \\ & \\ & 2 \sqrt{b} \cos{(\omega) A e^{\alpha(k+1)}} \cdot \cos{(\omega (k+1)+p)} \\ & - b \cdot A \cdot e^{\alpha \cdot k} \cdot \cos{(\omega k + p)} \\ \end{align}

Division by $A \cdot e^{\alpha \cdot k}$ and application of addition theorems for cos-functions gives

\begin{align} & e^{2 \alpha} [\cos{(\omega k + p)} \cos{(2 \omega)}-\sin{(\omega k + p)} \sin{(2 \omega)} ] = \\ \\ & 2 e^{\alpha}\sqrt{b}\cos{\omega} [\cos{(\omega k + p)} \cos{(\omega)}-\sin{(\omega k + p)} \sin{(\omega)}] \\ & - b \cos{(\omega k+p)} \\ \end{align}

Rearranging terms

\begin{align} 0=& & \cos{(\omega k + p)} [e^{2 \cdot \alpha} \cos{(2 \omega)}-2e^\alpha \sqrt b \cos ^2 {\omega} +b ]- \\ & & \sin{(\omega k + p)} [e^{2 \cdot \alpha} \sin{(2 \omega)}-2e^\alpha \sqrt b \sin {\omega}\sin {\omega} ]= \\ \\ 0= & & \cos{(\omega k + p)} [e^{2 \cdot \alpha} \cos{(2 \omega)}-2e^\alpha \sqrt b \cos ^2 {\omega} +b ]- \\ & & \sin{(\omega k + p)}\sin{(2 \omega)} \cdot [e^{2 \alpha } -e^{ \alpha }\sqrt b ] \\ \end{align}

The last equation can only hold for any k if the expressions in the squared brackets equal 0:

$$0= e^{2 \alpha } -e^{ \alpha }\sqrt b$$
$$\Rightarrow e^{\alpha } = \sqrt b$$

[ - ]
Comment by July 2, 2015 That is how technical articles should be, straight into core subject, no messing about.

Kaz
[ - ]
Comment by July 3, 2015 What has ? got to do with oversampling as written in the sentence "For high oversampling rates, i.e. small resulting ?". Intuitively oversampling should give better results shouldn't it? Why is the fitting best for "near half the Nyquist rate"?
[ - ]
Comment by July 4, 2015 Hi jin_,

$\omega$ is a normalized frequency, it is related to the sample frequency. So for a given signal, if the sample rate gets higher, the $\omega$ is going down. From the linear system of equations the $\cos \omega$ term is calculated. This function is flat for small $\omega$, a small change in the $\cos \omega$ will result in a big change for $\omega$. The situation ist best near $\pi / 2$, because the $\cos \omega$ function ist steepest there, and this sample rate is corresponding to Nyquist / 2 . And yes, higher oversampling rate gives better results, you use all samples, but you have to rearrange the samples. Take a look at the Matlab-script.

Cheers
Detlef
[ - ]
Comment by September 21, 2015 Hello Dr. Amberg. In my efforts to understand your blog I searched the MathWorks Exchange web page for 'matlab fit damped sine wave'. In doing so I received the "Your search returned No results" message. Just so you know, searching on 'Fit a damped sine wave' worked fine.

Regards,
[-Rick-]
[ - ]
Comment by September 22, 2015 Hi,
thanks for the hint. For Google 'matlab fit damped sine wave', the MathWork Exchange link is the first hit, at least for my Google results. I hope, that your efforts to understand my blog were crowned with success. I'd be more than glad to give some additional clarifications.

Cheers
Detlef

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.