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  
$$

Sign in