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
$$
- Comments
- Write a Comment Select to add a comment

Kaz


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

Regards,
[-Rick-]

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.
Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.