#LyX 2.3 created this file. For more info see http://www.lyx.org/ \lyxformat 544 \begin_document \begin_header \save_transient_properties true \origin unavailable \textclass memoir \begin_preamble \input{tex/preamble_article.tex} \end_preamble \options twoside,final \use_default_options true \maintain_unincluded_children false \language american \language_package babel \inputencoding utf8 \fontencoding global \font_roman "libertine" "Linux Libertine O" \font_sans "default" "default" \font_typewriter "default" "default" \font_math "libertine-ntxm" "auto" \font_default_family default \use_non_tex_fonts true \font_sc false \font_osf false \font_sf_scale 100 100 \font_tt_scale 100 100 \use_microtype false \use_dash_ligatures false \graphics default \default_output_format default \output_sync 1 \output_sync_macro "\synctex=1" \bibtex_command default \index_command default \paperfontsize 10 \spacing single \use_hyperref true \pdf_author "Dr.ir. J.A. de Jong - ASCEE" \pdf_bookmarks true \pdf_bookmarksnumbered false \pdf_bookmarksopen false \pdf_bookmarksopenlevel 1 \pdf_breaklinks true \pdf_pdfborder true \pdf_colorlinks true \pdf_backref false \pdf_pdfusetitle true \papersize a4paper \use_geometry true \use_package amsmath 1 \use_package amssymb 1 \use_package cancel 1 \use_package esint 1 \use_package mathdots 1 \use_package mathtools 1 \use_package mhchem 1 \use_package stackrel 1 \use_package stmaryrd 1 \use_package undertilde 1 \cite_engine biblatex \cite_engine_type authoryear \biblio_style plain \biblatex_bibstyle authoryear \biblatex_citestyle numeric \use_bibtopic false \use_indices false \paperorientation portrait \suppress_date false \justification true \use_refstyle 0 \use_minted 0 \index Index \shortcut idx \color #008000 \end_index \leftmargin 3cm \topmargin 3cm \rightmargin 2.5cm \bottommargin 3.5cm \headsep 1cm \secnumdepth 3 \tocdepth 3 \paragraph_separation skip \defskip smallskip \is_math_indent 0 \math_numbering_side default \quotes_style english \dynamic_quotes 0 \papercolumns 1 \papersides 1 \paperpagestyle default \tracking_changes false \output_changes false \html_math_output 0 \html_css_as_file 0 \html_be_strict false \end_header \begin_body \begin_layout Standard \begin_inset ERT status open \begin_layout Plain Layout \backslash date{ \backslash today} \end_layout \begin_layout Plain Layout % \backslash begin{center} \end_layout \begin_layout Plain Layout % \backslash includegraphics{/home/anne/nextcloud/template_huisstijl/lyx/%ascee_beeldmerk_wit hacr.eps} \end_layout \begin_layout Plain Layout % \backslash end{center} \end_layout \end_inset \end_layout \begin_layout Title LASP \begin_inset Newline newline \end_inset \size normal Library for Acoustic Signal Processing \end_layout \begin_layout Author \series bold J.A. de Jong \begin_inset Formula $^{1}$ \end_inset \end_layout \begin_layout Standard \begin_inset VSpace medskip \end_inset \end_layout \begin_layout Standard \align center \begin_inset Graphics filename /home/anne/nextcloud/templates_huisstijl/lyx/ascee_beeldmerk_withacr.eps width 45text% \end_inset \end_layout \begin_layout Standard \align center \size small \begin_inset Formula $^{1}$ \end_inset ASCEE, Máximastraat 1, 7442 NW Nijverdal, info@ascee.nl \end_layout \begin_layout Standard \begin_inset VSpace medskip \end_inset \end_layout \begin_layout Standard \align left \begin_inset Tabular \begin_inset Text \begin_layout Plain Layout \begin_inset ERT status collapsed \begin_layout Plain Layout \backslash arrayrulecolor{asceelightblue} \end_layout \begin_layout Plain Layout \backslash midrule[2pt] \end_layout \end_inset Internal document ID: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 31416 \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout External document ID: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Document status: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Draft \end_layout \end_inset \end_inset \end_layout \begin_layout Standard \align center \size small \begin_inset VSpace medskip \end_inset \end_layout \begin_layout Standard \begin_inset Note Note status open \begin_layout Plain Layout \align left \begin_inset Tabular \begin_inset Text \begin_layout Plain Layout \begin_inset ERT status collapsed \begin_layout Plain Layout \backslash arrayrulecolor{asceelightblue} \end_layout \begin_layout Plain Layout \backslash midrule[2pt] \end_layout \end_inset Internal document ID: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout External document ID: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Document status: \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Draft \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset Note Note status open \begin_layout Plain Layout \lang dutch Revisiehistorie: \end_layout \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset Note Note status open \begin_layout Plain Layout \lang dutch \begin_inset ERT status open \begin_layout Plain Layout \backslash today \end_layout \end_inset : rev. 1 \begin_inset Newline newline \end_inset \begin_inset ERT status open \begin_layout Plain Layout \backslash fourteendaysahead \end_layout \end_inset : rev. 2 \end_layout \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset ERT status collapsed \begin_layout Plain Layout \backslash \backslash \end_layout \begin_layout Plain Layout \backslash bottomrule[2pt]% \end_layout \end_inset \end_layout \end_inset \end_inset \end_layout \end_inset \end_layout \begin_layout Standard \begin_inset CommandInset toc LatexCommand tableofcontents \end_inset \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Plain Layout \backslash thispagestyle{empty} \end_layout \begin_layout Plain Layout % Optionally: set this document to confidential \end_layout \begin_layout Plain Layout % \backslash confidential \end_layout \end_inset \end_layout \begin_layout Chapter Power spectra estimation \end_layout \begin_layout Section Fourier transform vs discrete Fourier transform \end_layout \begin_layout Standard Our definition of the Fourier transform for the real-valued function \begin_inset Formula $x(t)$ \end_inset with unit \begin_inset Formula $U$ \end_inset is: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} X(f)=\int\limits _{t=-\infty}^{\infty}x(t)\exp\left(-i2\pi ft\right)\mathrm{d}t, \end{equation} \end_inset such that: \begin_inset Formula \begin{equation} x(t)=\int\limits _{f=-\infty}^{\infty}X(f)\exp\left(2\pi ift\right)\mathrm{d}f \end{equation} \end_inset \end_layout \begin_layout Standard Using the radian frequency \begin_inset Formula $\omega$ \end_inset , a scaling factor should be used in front: \begin_inset Formula \begin{equation} x(t)=\frac{1}{2\pi}\int\limits _{\omega=-\infty}^{\infty}X(\omega)\exp\left(i\omega t\right)\mathrm{d}\omega, \end{equation} \end_inset where \begin_inset Formula \begin{equation} X(\omega)=\int\limits _{\omega=-\infty}^{\infty}x(t)\exp\left(-i\omega t\right)\mathrm{d}t, \end{equation} \end_inset \end_layout \begin_layout Standard Our definition of the \bar under power spectral density \bar default is: \begin_inset Formula \begin{equation} P_{x}=\underbrace{\lim_{T\to\infty}\frac{1}{T}\int\limits _{t=-T}^{T}x^{2}(t)\mathrm{d}t}_{\mathrm{Signal\,power}}=E\left[x^{2}(t)\right]\equiv\int\limits _{f=-\infty}^{\infty}S_{xx}(f)\mathrm{d}f \end{equation} \end_inset \end_layout \begin_layout Standard From Parseval's theorem, we know that \begin_inset Formula \begin{equation} S_{xx}(f)=\lim_{T\to\infty}\frac{1}{T}X(f)X^{*}(f). \end{equation} \end_inset Hence for signals for which the Fourier transform formally exist, the power spectral density is zero. In practice, signal \end_layout \begin_layout Standard \begin_inset Note Note status collapsed \begin_layout Plain Layout Filling in: \end_layout \begin_layout Plain Layout \begin_inset Formula \[ \lim_{T\to\infty}\frac{1}{T}\int\limits _{t=-T}^{T}\left[\int\limits _{\omega=-\infty}^{\infty}X(\omega)\exp\left(i\omega t\right)\mathrm{d}\omega\int\limits _{\omega=-\infty}^{\infty}X^{*}(\omega)\exp\left(-i\omega t\right)\mathrm{d}\omega\right]\mathrm{d}t=\int\limits _{\omega=-\infty}^{\infty}S_{xx}(\omega)\mathrm{d}\omega \] \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula \[ \lim_{T\to\infty}\frac{1}{T}\int\limits _{t=-T}^{T}\left[\int\limits _{\omega=-\infty}^{\infty}X(\omega)\exp\left(i\omega t\right)\mathrm{d}\omega\int\limits _{\omega=-\infty}^{\infty}X^{*}(\omega)\exp\left(-i\omega t\right)\mathrm{d}\omega\right]\mathrm{d}t=\int\limits _{\omega=-\infty}^{\infty}S_{xx}(\omega)\mathrm{d}\omega \] \end_inset \end_layout \end_inset \begin_inset Formula \[ \] \end_inset \end_layout \begin_layout Standard Results in: \end_layout \begin_layout Standard Plancheler theorem: \begin_inset Formula \begin{equation} \int\limits _{t=-\infty}^{\infty}x(t)^{2}\mathrm{d}t=\int\limits _{\omega=-\infty}^{\infty}X(\omega)X^{*}(\omega)\mathrm{d}\omega \end{equation} \end_inset \end_layout \begin_layout Section Estimation of power spectra \end_layout \begin_layout Itemize Sample frequency: \begin_inset Formula $f_{s}=1/\Delta t$ \end_inset \end_layout \begin_layout Itemize \begin_inset Formula $N_{\mathrm{DFT}}$ \end_inset The number of samples taking into the FFT (nfft) \end_layout \begin_layout Standard Frequency resolution: \begin_inset Formula \begin{equation} \Delta f=\frac{f_{s}}{N_{\mathrm{DFT}}}, \end{equation} \end_inset i.e. the smallest frequency that fits into the measured time \begin_inset Formula $T=N_{\mathrm{DFT}}\Delta t$ \end_inset . \end_layout \begin_layout Standard Our definition of the DFT: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} X[k]=\mathrm{DFT}\left(x_{n}\right)=\sum_{n=0}^{N_{\mathrm{DFT}}-1}x_{n}e^{-2\pi ikn/N_{\mathrm{DFT}}},\label{eq:dft_definition} \end{equation} \end_inset such that \begin_inset Formula \begin{equation} x[n]=\mathrm{iDFT}\left(X[k]\right)=\frac{1}{N_{\mathrm{DFT}}}\sum_{k=0}^{N_{\mathrm{DFT}}-1}X[k]e^{2\pi ikn/N_{\mathrm{DFT}}}. \end{equation} \end_inset The advantage of this definition is that it preserves the duality between transfer function and impulse response. However, for proper scaling of signal (power) spectra, it requires more inspection. \end_layout \begin_layout Standard Using the definition of Eq. \begin_inset space ~ \end_inset \begin_inset CommandInset ref LatexCommand ref reference "eq:dft_definition" \end_inset , using the real-valued DFT (DFTR), we directly obtain the positive frequency half-spectrum. \end_layout \begin_layout Standard Parseval's theorem states: \begin_inset Formula \begin{equation} \sum_{n=0}^{N_{\mathrm{DFT}}-1}x[n]y[n]^{*}=\frac{1}{N}\sum_{k=0}^{N_{\mathrm{DFT}}-1}X[k]Y[k]^{*} \end{equation} \end_inset \end_layout \begin_layout Standard The average signal power is defined as \begin_inset Formula \begin{equation} P(x[n])=\frac{1}{N_{\mathrm{DFT}}}\sum_{n=0}^{N_{\mathrm{DFT}}-1}x[n]^{2} \end{equation} \end_inset which can be written in frequency domain, using Parseval's theorem as: \begin_inset Formula \begin{equation} P(x[n])=\frac{1}{N_{\mathrm{DFT}}^{2}}\sum_{k=0}^{N_{\mathrm{DFT}}-1}X[k]X[k]^{*}. \end{equation} \end_inset Hence the signal power as a function of frequency is \begin_inset Formula \begin{equation} P_{k}=\frac{1}{N_{\mathrm{DFT}}^{2}}\left\Vert X[k]\right\Vert ^{2}, \end{equation} \end_inset such that \begin_inset Formula $P_{k}$ \end_inset is the signal power in frequency bin \begin_inset Formula $k$ \end_inset . Now, the power spectral density is defined as the signal power per unit frequency, for which the signal power needs to be divided by the frequency resolution. Hence \begin_inset Formula \begin{equation} S_{x}=\frac{1}{N_{\mathrm{DFT}}f_{s}}\left\Vert X[k]\right\Vert ^{2} \end{equation} \end_inset \end_layout \begin_layout Subsection Single sided amplitude spectra \end_layout \begin_layout Standard If \begin_inset Formula $X_{s}$ \end_inset is defined such that \begin_inset Formula \begin{equation} x[n]=\frac{1}{\sqrt{2}}X_{s}[0]+\sum_{k=1}^{N_{\mathrm{DFT}}/2-1}X_{s}[k]e^{2\pi ikn/N_{\mathrm{DFT}}}+\frac{1}{\sqrt{2}}X_{s}[N_{\mathrm{DFT}}/2]e^{i\pi n}, \end{equation} \end_inset then the value of \begin_inset Formula $X_{s}$ \end_inset directly corresponds to the amplitude of the sinusoid, for all frequencies not equal to the DC and Nyquist rate. \begin_inset Formula $N_{\mathrm{DFT}}/2$ \end_inset is an integer division which rounds down to the nearest integer. From this definition, we can conclude that \begin_inset Formula \begin{eqnarray} X_{s}[0] & =\frac{\sqrt{2}}{N_{\mathrm{DFT}}}X[0] & \mathrm{for}\,\,\,k=0\\ X_{s}[k] & =\frac{2}{N_{\mathrm{DFT}}}X[k] & \mathrm{for}\,\,\,0 \begin_inset Text \begin_layout Plain Layout Dimension \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Value \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Lower frequency \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout Upper frequency \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset Formula $s_{1}$ \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 430 \begin_inset space ~ \end_inset mm \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 40 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 320 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset Formula $s_{2}$ \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 190 \begin_inset space ~ \end_inset mm \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 90 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 724 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout \begin_inset Formula $s_{3}$ \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 34 \begin_inset space ~ \end_inset mm \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 491 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \begin_inset Text \begin_layout Plain Layout 4047 \begin_inset space ~ \end_inset Hz \end_layout \end_inset \end_inset \end_layout \begin_layout Subsection Microphone switching technique for estimating the transfer function corrections \end_layout \begin_layout Itemize Normal measurement, mic 0 at position A, mic 1 at position B \end_layout \begin_layout Itemize Switched measurement, mic 0 at position B, mic 1 at position A \end_layout \begin_layout Standard Definitions: \end_layout \begin_layout Itemize \begin_inset Formula $K_{0}$ \end_inset : Microphone calibration correction factor for mic 0. Such that \begin_inset Formula $p_{0}=K_{0}\tilde{p}_{0}$ \end_inset , where \begin_inset Formula $\hat{p}_{0}$ \end_inset is the measured microphone pressure, and \begin_inset Formula $p_{0}$ \end_inset the actual pressure at the measurement position. \end_layout \begin_layout Itemize \begin_inset Formula $K_{1}$ \end_inset : Microphone calibration correction factor for mic 1. Such that \begin_inset Formula $p_{1}=K_{1}\tilde{p}_{1}$ \end_inset , where \begin_inset Formula $\hat{p}_{1}$ \end_inset is the measured microphone pressure, and \begin_inset Formula $p_{0}$ \end_inset the actual pressure at the measurement position. \end_layout \begin_layout Standard We are only able to measure and estimate cross-spectrum (or equivalently the cross-spectral density): \begin_inset Formula \begin{equation} \tilde{C}_{ij}=\tilde{p}_{i}\tilde{p}_{j}^{*}, \end{equation} \end_inset from which the transfer functions can be estimated. We require to correct these transfer functions for the relative microphone calibration. The final quantity of interest is often the acoustic pressure transfer function for the two mic's in the impedance tube: \begin_inset Formula \begin{equation} G_{AB}=\frac{p_{B}}{p_{A}} \end{equation} \end_inset \end_layout \begin_layout Standard he microphone switching method is used to estimate the calibration constant that should be used to estimate the transfer function from the measured CPS's. In measurement 1: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} G_{AB}^{(1)}=\frac{p_{1}^{(1)}}{p_{0}^{(1)}}\approx\frac{K_{1}}{K_{0}}\frac{\tilde{C}_{10}^{(1)}}{\tilde{C}_{00}^{(1)}}, \end{equation} \end_inset From the second measurement in switched configuration, the estimation of \begin_inset Formula $G_{AB}$ \end_inset yields: \begin_inset Formula \begin{equation} G_{AB}^{(2)}=\frac{p_{0}^{(2)}}{p_{1}^{(2)}}\approx\frac{K_{0}}{K_{1}}\frac{\tilde{C}_{01}^{(2)}}{\tilde{C}_{11}^{(2)}}, \end{equation} \end_inset \end_layout \begin_layout Standard Equating both expressions yields an exact expression for the calibration correction factor: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $G_{AB}^{(1)}=G_{AB}^{(2)}$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\left(\frac{K_{1}}{K_{0}}\right)^{2}=\frac{\tilde{C}_{00}^{(1)}}{\tilde{C}_{10}^{(1)}}\frac{\tilde{C}_{01}^{(2)}}{\tilde{C}_{11}^{(2)}}$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\frac{K_{1}}{K_{0}}=\sqrt{\frac{\tilde{C}_{01}^{(2)}}{\tilde{C}_{11}^{(2)}}\frac{\tilde{C}_{00}^{(1)}}{\tilde{C}_{10}^{(1)}}}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} K\equiv\frac{K_{1}}{K_{0}}=\sqrt{\frac{\tilde{C}_{01}^{(2)}}{\tilde{C}_{11}^{(2)}}\frac{\tilde{C}_{00}^{(1)}}{\tilde{C}_{10}^{(1)}}}. \end{equation} \end_inset Then, for all measurements beyond the calibration, we can write for \begin_inset Formula $G_{AB}$ \end_inset : \begin_inset Formula \begin{equation} G_{AB}=\frac{p_{B}}{p_{A}}=KG_{01}, \end{equation} \end_inset where \begin_inset Formula $G_{01}$ \end_inset is the measured transfer function from mic 0 to mic 1 ( \begin_inset Formula $G_{01}=p_{1}/p_{0}$ \end_inset ). \end_layout \begin_layout Section Method for computing the sample impedance, absorption and reflection coefficient \end_layout \begin_layout Standard \begin_inset Float figure wide false sideways false status open \begin_layout Plain Layout \noindent \align center \begin_inset Graphics filename img/imptube_meas_setups.pdf width 100text% \end_inset \end_layout \begin_layout Plain Layout \begin_inset Caption Standard \begin_layout Plain Layout Configurations for acoustic sample testing \end_layout \end_inset \begin_inset CommandInset label LatexCommand label name "fig:sampletesting_configs" \end_inset \end_layout \end_inset \end_layout \begin_layout Standard Fig. \begin_inset space ~ \end_inset \begin_inset CommandInset ref LatexCommand ref reference "fig:sampletesting_configs" \end_inset shows two configurations that can be used for acoustic sample testing. The top configuration is used for relatively open samples. \end_layout \begin_layout Subsubsection Open samples - configuration 1 \end_layout \begin_layout Standard For both microphones on the same side of the sample, the reflection coefficient at the position of microphone \begin_inset Formula $A$ \end_inset can be evaluated as: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} R_{A}\equiv R(x=0)=\frac{G_{AB}-e^{-iks}}{e^{iks}-G_{AB}}, \end{equation} \end_inset where \begin_inset Formula $k$ \end_inset is the wave number, and \begin_inset Formula $s$ \end_inset is the microphone spacing. The reflection coefficient rotates in phase going to the position of the sample: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} R=R_{A}\exp\left(2ik\left(s+d_{1}\right)\right). \end{equation} \end_inset The absorption coefficient \begin_inset Formula $\alpha$ \end_inset is: \begin_inset Formula \begin{equation} \alpha=1-\left|R\right|^{2} \end{equation} \end_inset \end_layout \begin_layout Standard The impedance at the sample location ( \begin_inset Formula $x=s+d_{1}$ \end_inset ) can be computed as: \begin_inset Formula \begin{equation} z=z_{0}\frac{1+R}{1-R} \end{equation} \end_inset \end_layout \begin_layout Standard For an arbitrary back side impedance, this impedance is due to the sample, and the back side impedance. \end_layout \begin_layout Standard For a thin sample (w.r.t. the wavelength) with back cavity, the sample impedance can be computed, by knowing the distance behind the sample ( \begin_inset Formula $d_{2}$ \end_inset ). We know for a fixed back cavity that, we set: \begin_inset Formula \begin{equation} \Delta p_{s}=Z_{s}U, \end{equation} \end_inset \end_layout \begin_layout Standard The impedance on the back side of the sample, for a closed cavity can be computed by assuming 100% reflection at the end of the back cavity: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $p=\cos\left(k\left(L-x\right)\right)$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $i\omega\rho_{0}U=-\frac{\partial p}{\partial x}\Rightarrow U=\frac{ik}{\omega\rho_{0}}\sin\left(k\left(L-x\right)\right)$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula \[ Z|_{x=0}=\frac{p}{U}|_{x=0}=\frac{\cos\left(kL\right)}{\frac{ik}{\omega\rho_{0}}\sin\left(kL\right)}=-iz_{0}\cot\left(kL\right) \] \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $R(x)=R|_{x=0}\exp\left(2ikx\right)$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $R(-L)=R_{0}\exp\left(-2ikL\right)$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $R_{0}=1$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $z=z_{0}\frac{1+R}{1-R}=z_{0}\frac{1+\exp\left(-2ikL\right)}{1-\exp\left(-2ikL\right)}=z_{0}\frac{\exp\left(ikL\right)+\exp\left(-ikL\right)}{\exp\left(ikL\right)-\exp\left(-ikL\right)}=-z_{0}i\frac{\cos\left(kL\right)}{\sin\left(kL\right)}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} z_{c}=-iz_{0}\cot\left(kd_{2}\right), \end{equation} \end_inset Using that, we are able to compute the sample impedance (jump impedance) as: \begin_inset Formula \begin{equation} z_{s}=2z_{0}\frac{1-Re^{2ikd_{2}}}{\left(R-1\right)\left(e^{2ikd_{2}}-1\right)} \end{equation} \end_inset \end_layout \begin_layout Subsection Error analysis \end_layout \begin_layout Standard Types of error: \end_layout \begin_layout Itemize Positioning errors, both sample and microphones \end_layout \begin_deeper \begin_layout Itemize Bias and random \end_layout \end_deeper \begin_layout Itemize Model errors (speed of sound, transport parameters) \end_layout \begin_layout Itemize Sensor errors \end_layout \begin_layout Itemize Nonlinearities \end_layout \begin_layout Itemize Noise \end_layout \begin_layout Chapter Digital signal processing \end_layout \begin_layout Section Filter bank design \end_layout \begin_layout Subsection FIR Filter \end_layout \begin_layout Standard An FIR filter performs the operation: \begin_inset Formula \begin{equation} y[n]=\boldsymbol{h}*\boldsymbol{x}[n]\equiv\sum_{m=0}^{N-1}h[m]x[n-m]. \end{equation} \end_inset This operation can be implemented in time domain, or in frequency domain. Notably the \begin_inset Quotes eld \end_inset fast convolution form \begin_inset Quotes erd \end_inset is of interest, as it requires fewer operations for high \begin_inset Formula $N$ \end_inset (smaller complexity). \end_layout \begin_layout Subsubsection Fast convolution: overlap-save method \end_layout \begin_layout Standard Definitions: \end_layout \begin_layout Itemize FIR Filter length: \begin_inset Formula $N$ \end_inset \end_layout \begin_layout Itemize Block length of inputs: \begin_inset Formula $L$ \end_inset \end_layout \begin_layout Standard Limitations: \end_layout \begin_layout Itemize Time sample block size \begin_inset Formula $L$ \end_inset should be larger than the filter order. \end_layout \begin_layout Standard \begin_inset Float figure wide false sideways false status open \begin_layout Plain Layout \align center \begin_inset Graphics filename img/overlap_save.pdf width 70text% \end_inset \begin_inset Caption Standard \begin_layout Plain Layout Overlap-save algorithm \end_layout \end_inset \begin_inset CommandInset label LatexCommand label name "fig:schematic_overlapsave" \end_inset \end_layout \end_inset \end_layout \begin_layout Standard The yellow blocks are saved for the next block. \end_layout \begin_layout Section Digital filters \end_layout \begin_layout Subsection \begin_inset Formula $z$ \end_inset -transform \end_layout \begin_layout Standard \begin_inset Formula $z$ \end_inset -transform of a sequence \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} Z\left(h[n]\right)=\sum_{n=-\infty}^{\infty}h[n]z^{-n} \end{equation} \end_inset \begin_inset Formula $z$ \end_inset -transform of the discrete impulse \begin_inset Formula \begin{equation} Z\left(\delta[n]\right)=z^{0}=1 \end{equation} \end_inset \begin_inset Formula $z$ \end_inset -transform of the unit step function \begin_inset Formula \begin{equation} Z\left(u[z]\right)=\frac{z}{z-1} \end{equation} \end_inset \begin_inset Formula $z$ \end_inset -transform of the exponentially decaying function \begin_inset Formula \begin{equation} f(t)=e^{-at}, \end{equation} \end_inset is discrete \begin_inset Formula \begin{equation} f[n]=e^{-anT}, \end{equation} \end_inset where \begin_inset Formula $T=f_{s}^{-1}$ \end_inset . \begin_inset Formula \begin{equation} Z\left(f[n]\right)=\frac{z}{z-e^{-aT}} \end{equation} \end_inset \end_layout \begin_layout Subsection DSP Equation \end_layout \begin_layout Standard Sample time is \begin_inset Formula $nT$ \end_inset , where \begin_inset Formula $T$ \end_inset is the sampling period (inverse sampling frequency). \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} y[n]=\sum_{m=1}^{M}a_{m}y[n-m]+\sum_{p=0}^{P}b_{p}x[n-p] \end{equation} \end_inset In the \begin_inset Formula $z$ \end_inset -domain, this equation can be written as \begin_inset Formula \begin{equation} Y[z]=\frac{\sum\limits _{p=0}^{P}b_{p}z^{-p}}{1-\sum\limits _{m=1}^{M}a_{m}z^{-m}}X[z]=H[z]\cdot X[z] \end{equation} \end_inset Analog input signal frequency: \begin_inset Formula $\omega$ \end_inset . Scaled frequency: \begin_inset Formula $\Omega=\omega T$ \end_inset . Frequency response of the scaled frequency goes from 0 to \begin_inset Formula $\pi$ \end_inset . Each DSP system has a frequency response that repeats at the sampling frequency ( \begin_inset Formula $\Omega=2\pi$ \end_inset ). The frequency response can be found from the \begin_inset Formula $z$ \end_inset -transform by filling in \begin_inset Formula $z=e^{sT}$ \end_inset , where \begin_inset Formula $s$ \end_inset is the Laplace variable. And filling in for \begin_inset Formula $s=i\omega$ \end_inset , to find: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} H\left(\omega\right)=H[z=e^{i\omega T}], \end{equation} \end_inset For the DSP equation: \begin_inset Formula \begin{equation} H(\omega)=\frac{\sum\limits _{p=0}^{P-1}a_{p}e^{-ip\Omega}}{1-\sum\limits _{m=1}^{M-1}b_{m}e^{-im\Omega}} \end{equation} \end_inset \end_layout \begin_layout Subsection FIR Filter design \end_layout \begin_layout Standard Transfer function of a FIR filter: \begin_inset Formula \begin{equation} T(e^{i\Omega})=\sum_{k=-N}^{N}a_{k}z^{-ik\Omega}\label{eq:fir_freq_response} \end{equation} \end_inset The left side of Eq. \begin_inset CommandInset ref LatexCommand ref reference "eq:fir_freq_response" \end_inset is the desired frequency response, The right side are the corresponding filter coefficients. Multiplying the LHS and RHS with \begin_inset Formula $e^{in\Omega}$ \end_inset and integrating from 0 to \begin_inset Formula $2\pi$ \end_inset , we find the following equation for the filter coefficient \begin_inset Formula $a_{n}$ \end_inset : \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $T(e^{i\Omega})=\sum_{k=-N}^{N}a_{k}z^{-i\Omega}$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\int_{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega=\int_{0}^{2\pi}\sum_{k=-N}^{N}a_{k}z^{-ik\Omega}e^{in\Omega}\mathrm{d}\Omega$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\int_{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega=\int_{0}^{2\pi}\sum_{k=-N}^{N}a_{k}z^{-ik\Omega}e^{in\Omega}\mathrm{d}\Omega$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\int_{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega=\int_{0}^{2\pi}\sum_{k=-N}^{N}a_{k}z^{i\Omega\left(n-k\right)}\mathrm{d}\Omega$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\int_{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega=\int_{0}^{2\pi}a_{n}\mathrm{d}\Omega$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\int_{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega=2\pi a_{n}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} a_{n}=\frac{1}{2\pi}\int\limits _{0}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega \end{equation} \end_inset This relation gives the coefficients of the \emph on noncausal \emph default form of the filter coefficients. As the coefficients are symmetrical, only for positive \begin_inset Formula $n$ \end_inset , the \begin_inset Formula $a_{n}$ \end_inset 's need to be computed. As the frequency spectrum is repeated and symmetrical around the Nyquist frequency, we can write this as: \begin_inset Formula \begin{equation} a_{n}=\frac{1}{2\pi}\left[\int\limits _{0}^{\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega+\int\limits _{\pi}^{2\pi}T(e^{i\Omega})e^{in\Omega}\mathrm{d}\Omega\right]. \end{equation} \end_inset If we specify a certain frequency response below the Nyquist frequency, and let the part above the Nyquist frequency be its mirror image, this can be written as \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $T\left(e^{i\left(\Omega+\pi\right)}\right)=T^{*}\left(e^{i\Omega}\right)$ \end_inset \end_layout \begin_layout Plain Layout Using \begin_inset Formula $k=\Omega-2\pi\Rightarrow\Omega=k+2\pi$ \end_inset , we can write the second integral as \end_layout \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{2\pi}\left[\int\limits _{0}^{\pi}Te^{in\Omega}\mathrm{d}\Omega+\int\limits _{k=-\pi}^{0}T(e^{i\left(k+2\pi\right)})e^{in\left(k+2\pi\right)}\mathrm{d}k\right].$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{2\pi}\left[\int\limits _{0}^{\pi}Te^{in\Omega}\mathrm{d}\Omega+\int\limits _{\Omega=-\pi}^{0}T(e^{i\left(\Omega+2\pi\right)})e^{in\left(\Omega+2\pi\right)}\mathrm{d}\Omega\right].$ \end_inset \end_layout \begin_layout Plain Layout Using the fact that \begin_inset Formula $e^{2in\pi}=1$ \end_inset for all integer \begin_inset Formula $n$ \end_inset : \end_layout \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{2\pi}\left[\int\limits _{-\pi}^{\pi}Te^{in\Omega}\mathrm{d}\Omega\right].$ \end_inset \end_layout \begin_layout Plain Layout But we know that \begin_inset Formula $T\left(e^{-i\Omega}\right)=T^{*}\left(e^{i\Omega}\right)$ \end_inset , such that: \end_layout \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{2\pi}\left[\int\limits _{0}^{\pi}T\left(e^{i\Omega}\right)e^{in\Omega}+T^{*}\left(e^{i\Omega}\right)e^{-in\Omega}\mathrm{d}\Omega\right]$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} a_{n}=\frac{1}{2\pi}\left[\int\limits _{0}^{\pi}T\left(e^{i\Omega}\right)e^{in\Omega}+T^{*}\left(e^{i\Omega}\right)e^{-in\Omega}\mathrm{d}\Omega\right] \end{equation} \end_inset \end_layout \begin_layout Subsubsection Ideal pass-band filter \end_layout \begin_layout Standard For an ideal band filter, this results in \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{2\pi}\int\limits _{\Omega_{l}}^{\Omega_{u}}e^{in\Omega}+e^{-in\Omega}\mathrm{d}\Omega=$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $a_{n}=\frac{1}{\pi}\int\limits _{\Omega_{l}}^{\Omega_{u}}\cos\left(n\Omega\right)\mathrm{d}\Omega=\frac{\sin\left(n\Omega\right)}{n\pi}|_{\Omega_{l}}^{\Omega_{u}}=\frac{\sin\left(n\Omega_{u}\right)-\sin\left(n\Omega_{l}\right)}{n\pi}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} a_{n}=\frac{1}{2\pi}\int\limits _{\Omega_{l}}^{\Omega_{u}}e^{in\Omega}+e^{-in\Omega}\mathrm{d}\Omega=\frac{\sin\left(n\Omega_{u}\right)-\sin\left(n\Omega_{l}\right)}{n\pi} \end{equation} \end_inset \end_layout \begin_layout Subsubsection Ideal low-pass filter \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} a_{n}=\frac{\sin\left(n\omega_{p}T\right)}{n\pi}, \end{equation} \end_inset where \begin_inset Formula $\omega_{p}$ \end_inset is the pass-band to stop-band transition frequency. \end_layout \begin_layout Subsubsection Gibbs phenomenon \end_layout \begin_layout Standard The abrupt change in filter gain from pass band to stop band results is pass-band ripple and . Leftover ripple at an abrupt transition in filter coefficients is about 9% of the gain amplitude. Windowing the filter coefficients: \end_layout \begin_layout Itemize Magnitude of the ripple decreases \end_layout \begin_layout Itemize Width of the transition band increases \end_layout \begin_layout Standard If the unwindowed filter coefficients are \begin_inset Formula $f[n]$ \end_inset , the windowed filter coefficients are \begin_inset Formula \begin{equation} a_{n}=w[n]f[n] \end{equation} \end_inset \end_layout \begin_layout Standard Hamming window is a bit better than Hann window. Hann window is an improvement over the Bartlett window. In terms of the transition band width. However, stop-band gain increases. \end_layout \begin_layout Section Band-pass digital filter design \end_layout \begin_layout Standard Low pass filter \end_layout \begin_layout Itemize \begin_inset Formula $g_{s,\mathrm{max}}$ \end_inset : maximum allowed gain in the stop-band \end_layout \begin_layout Itemize \begin_inset Formula $g_{p,\mathrm{max}}$ \end_inset : maximum allowed gain in the pass-band \end_layout \begin_layout Itemize \begin_inset Formula $g_{p,\mathrm{min}}$ \end_inset : minimum allowed gain in the pass-band \end_layout \begin_layout Subsection Stereo band pass filter display \end_layout \begin_layout Itemize Digital filter should be designed in frequency domain \end_layout \begin_layout Itemize The RMS of the output of each filter is computed over a period on once the center frequency of the band. \end_layout \begin_layout Standard Example: 1000 Hz octave band: pass-band is from 750 Hz to 1500 Hz. Use a Hamming window to reduce the pass and stop-band ripple. \end_layout \begin_layout Subsection First order digital high pass \end_layout \begin_layout Standard A simple digital high-pass filter can be implemented using: \begin_inset Formula \begin{equation} G(s)=\frac{\tau s}{1+\tau s}, \end{equation} \end_inset where \begin_inset Formula $\tau$ \end_inset is the \begin_inset Formula $-3$ \end_inset \begin_inset space ~ \end_inset dB time constant, as when \begin_inset Formula $\omega=\tau^{-1}$ \end_inset , \begin_inset Note Note status collapsed \begin_layout Plain Layout \end_layout \end_inset \begin_inset Formula $|G|=\frac{1}{\sqrt{2}}$ \end_inset . \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $2\left(\omega\tau\right)^{2}=1+\left(\omega\tau\right)^{2}$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\left(\omega\tau\right)^{2}=1$ \end_inset \end_layout \end_inset Such that the cut-on frequency is: \begin_inset Formula \begin{equation} \frac{1}{2\pi f_{c}}=\tau. \end{equation} \end_inset Filling in, we find: \begin_inset Formula \begin{equation} G(s)=\frac{\frac{s}{2\pi f_{c}}}{1+\frac{s}{2\pi f_{c}}}, \end{equation} \end_inset Applying the bilinear transform, we find: \begin_inset Note Note status open \begin_layout Plain Layout \begin_inset Formula $G(s)=\frac{2\pi f_{c}s}{1+2\pi f_{c}s},$ \end_inset \end_layout \begin_layout Plain Layout Bilinear: \begin_inset Formula $2f_{s}\frac{z-1}{z+1},$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\frac{2f_{s}\tau\left(z-1\right)}{2f_{s}\tau\left(z-1\right)+z+1}$ \end_inset \end_layout \begin_layout Plain Layout Readjust: \end_layout \begin_layout Plain Layout \begin_inset Formula $\frac{2f_{s}\tau\left(1-z^{-1}\right)}{2f_{s}\tau\left(1-z^{-1}\right)+1+z^{-1}}$ \end_inset \end_layout \begin_layout Plain Layout Make denominator a0 1: \end_layout \begin_layout Plain Layout \begin_inset Formula $\frac{\frac{2f_{s}\tau}{\left(1+2f_{s}\tau\right)}\left(1-z^{-1}\right)}{1+\frac{\left(1-2f_{s}\tau\right)}{\left(1+2f_{s}\tau\right)}z^{-1}}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} G[z]=\frac{\frac{2f_{s}\tau}{\left(1+2f_{s}\tau\right)}\left(1-z^{-1}\right)}{1+\frac{\left(1-2f_{s}\tau\right)}{\left(1+2f_{s}\tau\right)}z^{-1}} \end{equation} \end_inset \end_layout \begin_layout Section Sound level meter implementation \end_layout \begin_layout Subsection Time-weighted sound level \end_layout \begin_layout Standard Fast time-weighted sound level,of the A-weighted pressure signal \begin_inset Formula $p_{A}(t)$ \end_inset : \begin_inset Formula \begin{equation} L_{AF}=10\log_{10}\left(\frac{1}{\tau_{F}}\int\limits _{-\infty}^{t}p_{A}^{2}\left(\xi\right)e^{-\left(t-\xi\right)/\tau_{F}}\mathrm{d}\xi\right)-10\log_{10}\left(p_{\mathrm{ref}}^{2}\right), \end{equation} \end_inset where \begin_inset Formula $\tau_{F}$ \end_inset is the exponential time constant in seconds for the fast time weighting. Implementation suggestion: square the frequency-weighted input signal, and apply a single pole low-pass filter with one pole at \begin_inset Formula $-\tau_{F}^{-1}$ \end_inset . \end_layout \begin_layout Itemize Fast time weighting: \begin_inset Formula $\tau_{F}=0,125$ \end_inset s \end_layout \begin_layout Itemize Slow time weighting: \begin_inset Formula $\tau_{s}=1$ \end_inset s. \end_layout \begin_layout Itemize Impulse time weighting: \begin_inset Formula $\tau_{i}=35$ \end_inset \begin_inset space ~ \end_inset ms \end_layout \begin_layout Subsection Implementation of single pole low pass filter \end_layout \begin_layout Standard The time weighting is specified as: \begin_inset Formula \begin{equation} L_{\mathrm{AF}}(t)=10\log_{10}\left[\frac{1}{p_{\mathrm{ref}}^{2}}\frac{1}{\tau_{F}}\int_{-\infty}^{t}p_{\mathrm{A}}^{2}(t)\exp\left(-\frac{\left(t-\xi\right)}{\tau_{\mathrm{F}}}\right)\mathrm{d}\xi\right]. \end{equation} \end_inset The Laplace transform of the integral is: \begin_inset Formula \begin{equation} \mathcal{L}\left[\frac{1}{\tau_{F}}\int_{-\infty}^{t}p_{\mathrm{A}}^{2}\exp\left(-\frac{\left(t-\xi\right)}{\tau_{\mathrm{F}}}\right)\mathrm{d}\xi\right]=p_{\mathrm{A}}^{2}(s)\frac{1}{1+\tau_{F}s} \end{equation} \end_inset \end_layout \begin_layout Standard A single pole low pass filter has a frequency response of \begin_inset Formula \begin{equation} G_{\mathrm{splp}}=\frac{1}{1+\tau s}, \end{equation} \end_inset \end_layout \begin_layout Subsubsection Implementation using Bilinear transform \end_layout \begin_layout Standard we create a digital filter from this one using the bilinear transform: \begin_inset Formula \begin{equation} s\to2f_{s}\frac{z-1}{z+1}, \end{equation} \end_inset which yields the digital filter: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $G_{\mathrm{splp},d}=\frac{1}{1+\tau2f_{s}\frac{z-1}{z+1}}=\frac{z+1}{z+1+\tau2f_{s}\left(z-1\right)}$ \end_inset \end_layout \begin_layout Plain Layout Filling in \begin_inset Formula $z=\exp(i\omega/f_{s})$ \end_inset and setting \begin_inset Formula $\omega=0$ \end_inset \end_layout \begin_layout Plain Layout \series bold \begin_inset Formula $G_{\mathrm{splp},d}(\omega=0)=1$ \end_inset \end_layout \begin_layout Plain Layout ============== \end_layout \begin_layout Plain Layout \begin_inset Formula $G_{\mathrm{splp},d}=\frac{1+z^{-1}}{\left(1+2f_{s}\tau\right)+\left(1-2f_{s}\tau\right)z^{-1}}$ \end_inset \end_layout \begin_layout Plain Layout Normalizing such that \begin_inset Formula $a_{0}=1$ \end_inset : \end_layout \begin_layout Plain Layout \family roman \series medium \shape up \size normal \emph off \bar no \strikeout off \xout off \uuline off \uwave off \noun off \color none \begin_inset Formula $G_{\mathrm{splp},d}=\frac{\left(1+2\tau f_{s}\right)^{-1}\left(1+z^{-1}\right)}{1+\frac{\left(1-\tau2f_{s}\right)}{\left(1+2\tau f_{s}\right)}z^{-1}}$ \end_inset \end_layout \begin_layout Plain Layout – Check: filling in \begin_inset Formula $z=0$ \end_inset should give a unit gain: \end_layout \begin_layout Plain Layout \family roman \series medium \shape up \size normal \emph off \bar no \strikeout off \xout off \uuline off \uwave off \noun off \color none \begin_inset Formula $G_{\mathrm{splp},d}(z=0)=\frac{\left(1+2\tau f_{s}\right)^{-1}}{\frac{\left(1-\tau2f_{s}\right)}{\left(1+2\tau f_{s}\right)}}$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} G_{\mathrm{splp},d}=\frac{\left(1+2\tau f_{s}\right)^{-1}\left(1+z^{-1}\right)}{1+\frac{\left(1-\tau2f_{s}\right)}{\left(1+2\tau f_{s}\right)}z^{-1}}, \end{equation} \end_inset So its digital filter coefficients are: \begin_inset Formula \begin{align} \boldsymbol{b} & =\left[\begin{array}{ccc} \left(1+2\tau f_{s}\right)^{-1} & \left(1+2\tau f_{s}\right)^{-1} & 0\end{array}\right]^{\mathrm{T}}\\ \boldsymbol{a} & =\left[\begin{array}{ccc} 1 & \frac{\left(1-\tau2f_{s}\right)}{\left(1+2\tau f_{s}\right)} & 0\end{array}\right]^{\mathrm{T}} \end{align} \end_inset \end_layout \begin_layout Standard No correction for frequency warping has been done, as for all cases \begin_inset Formula $\tau f_{s}\gg1$ \end_inset . The output frequency of the sound level meter will be decimated to a sampling frequency, where the single pole low pass filter has a -20 dB point, such that aliasing is at max ~ 0.1 times any oscillation: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $\frac{1-0.01}{0.01}=\left|\tau s\right|\Rightarrow\left|s\right|=\frac{1}{\tau}\frac{1-0.01}{0.01}\Rightarrow f_{s,\mathrm{slm}}=$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} \left|G_{\mathrm{splp}}\right|=\left|\frac{1}{1+\tau s}\right|=-20\,\mathrm{dB}=0.01\Rightarrow s=\frac{1}{\tau}\frac{1-0.01}{0.01} \end{equation} \end_inset \end_layout \begin_layout Standard Hence the minus 20 \begin_inset space ~ \end_inset dB point lies at a sampling frequency of: \begin_inset Formula \begin{equation} f_{s,\mathrm{slm}}=\frac{1}{2\pi\tau}\frac{1-0.01}{0.01} \end{equation} \end_inset \end_layout \begin_layout Standard Then, the downsampling factor is set at: \begin_inset Formula \begin{equation} d=\left\lfloor \frac{f_{s}}{f_{s,\mathrm{slm}}}\right\rfloor , \end{equation} \end_inset where \begin_inset Formula $\left\lfloor \dots\right\rfloor $ \end_inset denotes the floor operation. \end_layout \begin_layout Standard Each subsample corresponds to \begin_inset Formula \begin{equation} n=o+id \end{equation} \end_inset \end_layout \begin_layout Standard Find the highest \begin_inset Formula $i$ \end_inset , that does not fit into \begin_inset Formula $N$ \end_inset anymore: \begin_inset Formula \[ id>N-o \] \end_inset \end_layout \begin_layout Subsubsection Implementation using the matched Z-transform method \begin_inset Formula \begin{equation} G_{\mathrm{splp},d}=\frac{1-\exp\left(-\left(f_{s}\tau\right)^{-1}\right)}{1-\exp\left(-\left(f_{s}\tau\right)^{-1}\right)z^{-1}} \end{equation} \end_inset \end_layout \begin_layout Standard By setting unity gain for \begin_inset Formula $\omega=0$ \end_inset , this results in the following simple difference equation: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula $Y[z]=\frac{1-\exp\left(\left(f_{s}\tau_{d}\right)^{-1}\right)}{1-\exp\left(\left(f_{s}\tau_{d}\right)^{-1}\right)z^{-1}}X[z]$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $y[n]=\exp\left(\left(f_{s}\tau_{d}\right)^{-1}\right)y[n-1]+\left[1-\exp\left(\left(f_{s}\tau_{d}\right)^{-1}\right)\right]x[n]$ \end_inset \end_layout \end_inset \begin_inset Formula \begin{align} y[n] & =\alpha y[n-1]+\left(1-\alpha\right)x[n],\\ \alpha & =\exp\left(-\frac{1}{f_{s}\tau}\right) \end{align} \end_inset \end_layout \begin_layout Standard This way: the number of dB's / s, that the level is decaying, after signal stop is: \begin_inset Note Note status collapsed \begin_layout Plain Layout \begin_inset Formula \[ y[n]=\alpha y[n-1] \] \end_inset \end_layout \begin_layout Plain Layout For \begin_inset Formula $\alpha>0$ \end_inset : \begin_inset Formula \[ |y[n]|=\alpha|y[n-1]| \] \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula \[ y[n]=y[0]\alpha^{n} \] \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula \[ L\left(y[n]\right)=20\log\left(|y[n]\right)=20\log\left(|y(0)|\right)+20\log\left(\alpha^{n}\right) \] \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula \begin{align*} L(y(t)) & =20\log\left(|y[n=\frac{t}{\Delta t}]\right)=20\log\left(|y(0)|\right)+\frac{t}{\Delta t}20\log\left(\alpha\right)\\ & =20\log\left(|y(0)|\right)+t\frac{1}{\Delta t}\underbrace{20\log\left(\alpha\right)} \end{align*} \end_inset \end_layout \end_inset \begin_inset Formula \begin{equation} \mathrm{Decay\,rate\,dB\,/\,s}\approx\frac{20}{\Delta t}\log\left(\alpha\right) \end{equation} \end_inset \begin_inset Note Note status collapsed \begin_layout Plain Layout And the other way around: \begin_inset Formula $20\log\left(\alpha\right)=d\Delta t$ \end_inset \end_layout \begin_layout Plain Layout \begin_inset Formula $\alpha=10^{\frac{d\Delta t}{20}}$ \end_inset \end_layout \end_inset \end_layout \begin_layout Subsection Time-averaged sound level \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} L_{A\mathrm{eq},T}=10\log_{10}\left(\frac{1}{p_{\mathrm{ref}}^{2}T}\int\limits _{t-T}^{t}p_{A}^{2}(\xi)\mathrm{d}\xi\right), \end{equation} \end_inset \end_layout \begin_layout Subsection C-weighting filter \end_layout \begin_layout Standard Linear amplitude scaling for the C-weighted frequency response: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} C(f)=C_{1000}^{-1}\left(\frac{f_{4}^{2}f^{2}}{\left(f^{2}+f_{1}^{2}\right)\left(f^{2}+f_{4}^{2}\right)}\right)^{2},\label{eq:C_weighting_freqresponse} \end{equation} \end_inset where \begin_inset Formula $C_{1000}$ \end_inset is the numerator of Eq. \begin_inset CommandInset ref LatexCommand ref reference "eq:C_weighting_freqresponse" \end_inset evaluated at 1 kHz. In this Eq. \begin_inset Formula \begin{equation} f_{1}=\sqrt{\left(\frac{-b-\sqrt{b^{2}-4c}}{2}\right)}\label{eq:f_1} \end{equation} \end_inset and \begin_inset Formula \begin{equation} f_{4}=\sqrt{\left(\frac{-b+\sqrt{b^{2}-4c}}{2}\right)}\label{eq:f_4} \end{equation} \end_inset , where \begin_inset Formula $c=f_{L}^{2}f_{H}^{2}$ \end_inset and \begin_inset Formula \begin{equation} b=\frac{1}{1-D}\left[f_{r}^{2}+\frac{f_{L}^{2}f_{H}^{2}}{f_{r}^{2}}-D\left(f_{L}^{2}+f_{H}^{2}\right)\right] \end{equation} \end_inset , with \begin_inset Formula $D=+\frac{1}{2}\sqrt{2}$ \end_inset . \begin_inset Formula $f_{r}=1$ \end_inset kHz, \begin_inset Formula $f_{L}=10^{1.5}$ \end_inset Hz and \begin_inset Formula $f_{H}=10^{3.9}$ \end_inset Hz. \end_layout \begin_layout Subsection A-weighting filter \end_layout \begin_layout Standard Linear amplitude scaling for the A-weighted frequency response: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} A(f)=A_{1000}^{-1}\frac{f_{4}^{2}f^{4}}{\left(f^{2}+f_{1}^{2}\right)\sqrt{\left(f^{2}+f_{2}^{2}\right)\left(f^{2}+f_{3}^{2}\right)}\left(f^{2}+f_{4}^{2}\right)},\label{eq:A_freq_norm} \end{equation} \end_inset in which \begin_inset Formula $f_{1}$ \end_inset and \begin_inset Formula $f_{4}$ \end_inset are defined in Eqs. \begin_inset CommandInset ref LatexCommand ref reference "eq:f_1" \end_inset and \begin_inset CommandInset ref LatexCommand ref reference "eq:f_4" \end_inset , respectively. And \begin_inset Formula \begin{equation} f_{2}=\frac{3-\sqrt{5}}{2}f_{A}, \end{equation} \end_inset and \begin_inset Formula \begin{equation} f_{3}=\frac{3+\sqrt{5}}{2}f_{A}, \end{equation} \end_inset where \begin_inset Formula $f_{A}=10^{2.45}$ \end_inset Hz. The transfer function of \begin_inset Formula $A$ \end_inset can be written as Equation \begin_inset CommandInset ref LatexCommand ref reference "eq:A_freq_norm" \end_inset can be rewritten to the following equivalent form \begin_inset CommandInset citation LatexCommand cite key "rimell_design_2015" literal "false" \end_inset : \begin_inset Formula \begin{equation} A(s)=K_{A}\frac{\omega_{4}^{2}s^{4}}{\left(s^{2}+\omega_{1}^{2}\right)\left(s+\omega_{2}\right)\left(s+\omega_{3}\right)\left(s^{2}+\omega_{4}^{2}\right)} \end{equation} \end_inset \end_layout \begin_layout Subsection Exact midband frequencies \end_layout \begin_layout Standard The midband frequencies of a (fractional) octave band filter are defined as: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} f_{x}=f_{r}\left[10^{\left(\frac{3}{10}\right)\left(\frac{x}{b}\right)}\right], \end{equation} \end_inset where \begin_inset Formula $f_{r}$ \end_inset is the reference frequency of 1000 Hz, \begin_inset Formula $10^{\left(3/10\right)}$ \end_inset is the nominal octave ratio for a base-10 system. \begin_inset Formula $b$ \end_inset is the step-width designator, \begin_inset Formula $b=3$ \end_inset for one-third-octave intervals. \end_layout \begin_layout Standard Octave ratio: Nominal frequency ratio of 2:1. Base ten system is preferred, where \begin_inset Formula \begin{equation} G_{10}=10^{3/10}\approx1.995 \end{equation} \end_inset , in base 2: \begin_inset Formula \begin{equation} G_{2}=2, \end{equation} \end_inset Bandwidth designator ( \begin_inset Formula $b$ \end_inset ). Exact midband frequencies: \begin_inset Formula \begin{align} f_{m} & =\left(G^{x/b}\right)f_{r},;\,b\,\mathrm{even}\\ f_{m} & =\left(G^{\left(2x+1\right)/(2b)}\right)f_{r},;\,b\,\mathrm{odd} \end{align} \end_inset where \begin_inset Formula $f_{r}$ \end_inset is the reference frequency (1 kHz). Bandedge frequencies: frequencies of the lower and upper edges of the passband of a bandpass filter suchthat the exact midband frequency is the geometric mean of the lower and upper bandedge frequencies. \begin_inset Formula \begin{align} f_{\ell} & =\left(G^{-1/\left(2b\right)}\right)f_{m}\\ f_{u} & =\left(G^{+1/\left(2b\right)}\right)f_{m} \end{align} \end_inset where \begin_inset Formula $f_{m}$ \end_inset is the exact midband frequency. \end_layout \begin_layout Subsection Nominal midband frequencies \end_layout \begin_layout Standard - See standard \end_layout \begin_layout Section Reverberation time \end_layout \begin_layout Standard Reverberation time is the time that is required to let the instantaneous sound pressure level drop with 60 db. Model: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} P^{2}(t)=P_{0}^{2}10^{-\alpha t} \end{equation} \end_inset Then at \begin_inset Formula $T_{60}$ \end_inset : \begin_inset Formula \begin{equation} \frac{P^{2}}{P_{0}^{2}}=10^{\frac{-60}{20}}=10^{-3} \end{equation} \end_inset \end_layout \begin_layout Standard Such that at \begin_inset Formula $T_{60}$ \end_inset , \begin_inset Formula $\alpha T_{60}=3$ \end_inset , hence \begin_inset Formula \begin{equation} T_{60}=\frac{3}{\alpha}. \end{equation} \end_inset \begin_inset Formula \begin{equation} P^{2}(t)=P_{0}^{2}10^{-\frac{3}{T_{60}}t} \end{equation} \end_inset \end_layout \begin_layout Standard Suppose we have the data \begin_inset Formula $\boldsymbol{P}^{2}$ \end_inset for given time instances \begin_inset Formula $\boldsymbol{t}$ \end_inset . Then we determine \begin_inset Formula $\alpha$ \end_inset and \begin_inset Formula $P_{0}^{2}$ \end_inset by minimizing: \begin_inset Formula \begin{equation} \left\Vert \boldsymbol{P}^{2}-P_{0}^{2}10^{-\alpha\boldsymbol{t}}\right\Vert ^{2}, \end{equation} \end_inset which is a nonlinear least squares problem \end_layout \begin_layout Standard Look at Schroeder-back integration! \end_layout \begin_layout Section Sweep signals \end_layout \begin_layout Standard A sweep signal is a signal which instantaneous frequency changes over time with a specific profile, the signal follows a sine function with a specific phase: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} c(t)=\sin\left(\phi(t)\right), \end{equation} \end_inset where \begin_inset Formula $\phi(t)$ \end_inset is the phase. For a sine wave with frequency \begin_inset Formula $f$ \end_inset , \begin_inset Formula $\phi(t)=2\pi ft$ \end_inset . The instantaneous frequency \begin_inset Formula $f(t)$ \end_inset is defined as \begin_inset Formula \begin{equation} f(t)=\frac{\mathrm{d}\phi}{\mathrm{d}t} \end{equation} \end_inset A linear sweep exhibits a profile where \begin_inset Formula $\frac{\mathrm{d}^{2}\phi}{\mathrm{d}t^{2}}=\mathrm{constant}$ \end_inset . For all the sweep implementations, we specify: \end_layout \begin_layout Itemize The lower frequency \begin_inset Formula $f_{L}$ \end_inset , \end_layout \begin_layout Itemize The upper frequency \begin_inset Formula $f_{U}$ \end_inset , \end_layout \begin_layout Itemize The sweep time of a single period \begin_inset Formula $T_{c}$ \end_inset , \end_layout \begin_layout Itemize The sweep profile, i.e. forward, backward or continuous \end_layout \begin_layout Itemize The sweep type (linear, exponential, hyperbolic) \end_layout \begin_layout Standard For the numerical implementation, the phase updating is done based on a forward Euler estimate of the derivative of the phase: \begin_inset Formula \begin{equation} \phi_{n+1}=\phi_{n}+2\pi f_{n}\Delta t,\qquad n=0...N-1 \end{equation} \end_inset where \begin_inset Formula $N=\left\lfloor T_{c}f_{s}\right\rfloor $ \end_inset , and \begin_inset Formula $\left\lfloor \dots\right\rfloor $ \end_inset denotes rounding down to the nearest integer. The input parameters are slightly adjusted, such that the sweep becomes periodic in the discrete time domain, with with \begin_inset Formula $C_{0}$ \end_inset -continuity, and a sign of the derivative pointing in the same direction. Note that for this type of \begin_inset Formula $C_{0}$ \end_inset -continuity, we require that \begin_inset Formula $\phi_{N}=2\pi K$ \end_inset , where \begin_inset Formula $K\in\mathbb{N}_{0}$ \end_inset (all natural numbers including 0). The boundary conditions are thus: \begin_inset Formula \begin{align} \phi_{0} & =0\\ \phi_{N} & =0\\ \frac{\mathrm{d}\phi_{0}}{\mathrm{d}n} & =\frac{\mathrm{d}\phi_{N}}{\mathrm{d}n} \end{align} \end_inset \end_layout \begin_layout Subsection Forward linear sweep \end_layout \begin_layout Standard For a forward linear sweep with period \begin_inset Formula $T_{c}$ \end_inset , the discrete sweep phase in a single period is defined as: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} f_{n}=f_{L}+\frac{n}{N}\left(f_{U}+\varepsilon-f_{L}\right).\label{eq:forward_linear} \end{equation} \end_inset This results in the following equation for the phase: \begin_inset Formula \begin{align} \phi_{n+1}=\phi_{n}+2\pi f_{L}\Delta t+2\pi\frac{n}{N}\Delta t\left(f_{U}+\varepsilon-f_{L}\right) \end{align} \end_inset Solving this recurrence relation yields: \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} \phi_{n}=2\pi\Delta t\left[f_{L}n+\frac{1}{2N}\left(n^{2}-n\right)\left(f_{U}+\varepsilon-f_{L}\right)\right] \end{equation} \end_inset \end_layout \begin_layout Standard We first make a proper estimate of \begin_inset Formula $K$ \end_inset , by setting \begin_inset Formula $\phi_{N}=2\pi K$ \end_inset : \begin_inset Formula \begin{equation} K=\left\lfloor \Delta t\left[f_{L}N+\frac{1}{2}\left(N-1\right)\left(f_{U}-f_{L}\right)\right]\right\rfloor , \end{equation} \end_inset which is used to set the correction \begin_inset Formula $\varepsilon$ \end_inset : \begin_inset Formula \begin{equation} \varepsilon=\frac{\frac{K}{\Delta t}-f_{L}N-\frac{1}{2}\left(N-1\right)\left(f_{U}-f_{L}\right)}{\frac{1}{2}\left(N-1\right)} \end{equation} \end_inset \end_layout \begin_layout Subsection Backward linear sweep \end_layout \begin_layout Standard The procedure for creating a backward linear sweep is similar to a forward linear sweep, only we replace \begin_inset Formula $f_{U}$ \end_inset with \begin_inset Formula $f_{L}$ \end_inset and vice versa. \end_layout \begin_layout Subsection Continuous linear sweep \end_layout \begin_layout Standard We define \begin_inset Formula $N_{f}$ \end_inset as \begin_inset Formula $\left\lfloor N/2\right\rfloor $ \end_inset , and \begin_inset Formula $N_{b}=N-N_{f}$ \end_inset . \end_layout \begin_layout Standard \begin_inset Formula \begin{equation} \phi_{n+1}=\begin{cases} \phi_{n}+2\pi\Delta t\left(f_{L}+\frac{n}{N_{f}}\left(f_{U}-f_{L}\right)\right) & 0\leq n\leq N_{f}\\ \phi_{n}+2\pi\Delta t\left(f_{U}+\varepsilon-\frac{n-N_{f}}{N_{b}}\left(f_{U}+\varepsilon-f_{L}\right)\right) & N_{f}