Next  |  Prev  |  Up  |  Top  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

Trapezoidal Rule for Numerical Integration

The velocity $ v(t)$ can be written as

\begin{eqnarray*}
v(t)=v(0)+\left(\int_{0}^{t}\dot v(\tau)d\tau\right)
\end{eqnarray*}

In particular,

\begin{eqnarray*}
v(nT) % &=& v(0) + \int_{0}^{nT}\vt(\tau)d\tau \\ [10pt]
&=& v(0) + \int_{0}^{(n-1)T}\dot v(\tau)d\tau + \int_{(n-1)T}^{nT}\dot v(\tau)d\tau \\ [10pt]
&=& v[(n-1)T]+\int_{(n-1)T}^{nT}\dot v(\tau)d\tau \\ [10pt]
&\approx& v[(n-1)T] +T\frac{\dot v[(n-1)T] + \dot v(nT)}{2}
\end{eqnarray*}

Bilinear Transform as Compensated BE/FE

In Newton's law $ f=m\dot v$ , look at the Backward Euler (BE) approximation of the time-derivative:

$\displaystyle f(t) \;=\;m\,\dot v\,\;\approx\;\,m\, \frac{v(t) - v(t-T)}{T}
$

We see there is a 1/2 sample delay in the first-order difference on the right. This misaligns the force $ f(t)$ and subsequent velocity by half a sample. A very simple delay compensation is to use a two-point average on the left:

$\displaystyle \frac{f(n)+f(n-1)}{2} \,\;\approx\;\,m\, \frac{v(n) - v(n-1)}{T}
$

The extra attenuation at high frequencies due to the two-point average actually helps. Taking the $ z$ transform:

$\displaystyle \frac{1+z^{-1}}{2}F(z) \,\;\approx\;\,m\, \frac{1-z^{-1}}{T} V(z)
$

or

$\displaystyle F(z) \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}\right) V(z)
$

which is the bilinear transform of $ F(s) = ms\,V(s)$ :

$\displaystyle \zbox{s \mapsto \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}
$

Frequency Warping is the Only Error

We have

$\displaystyle F(z) \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}\right) V(z)
$

using the bilinear transform (trapezoidal integration in the time domain)

Let's look along the unit circle in the $ z$ plane:

$\displaystyle \frac{F(e^{j\omega T})}{V(e^{j\omega T})} \,\;\approx\;\,m\, \left(\frac{2}{T} \frac{1-e^{-j\omega T}}{1+e^{-j\omega T}}\right)
\;=\;m\, j \left(\frac{2}{T} \tan\left(\frac{\omega T}{2}\right)\right)
$

Since the exact formula is $ F(e^{j\omega T})/V(e^{j\omega T}) = m\, j\omega$ , we can push all of the error into a frequency warping:

$\displaystyle \omega_d \mathrel{\stackrel{\mathrm{\Delta}}{=}}\frac{2}{T} \tan\left(\frac{\omega_a T}{2}\right)
$


Next  |  Prev  |  Up  |  Top  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

Download DigitizingNewton.pdf
Download DigitizingNewton_2up.pdf
Download DigitizingNewton_4up.pdf

``Introduction to Physical Signal Models'', by Julius O. Smith III, (From Lecture Overheads, Music 420).
Copyright © 2020-06-27 by Julius O. Smith III
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA  [Automatic-links disclaimer]