#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}