ASCEE
/
lasp_doc
4
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
lasp_doc/ASCEE_lasp.lyx

4839 lines
87 KiB
Plaintext
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#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
<lyxtabular version="3" rows="3" columns="2">
<features tabularvalignment="middle" tabularwidth="100text%">
<column alignment="left" valignment="top" width="0pt">
<column alignment="right" valignment="top" width="8cm">
<row>
<cell alignment="left" valignment="top" usebox="none">
\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
</cell>
<cell alignment="right" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
31416
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="left" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
External document ID:
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="left" valignment="top" bottomline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Document status:
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" bottomline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Draft
\end_layout
\end_inset
</cell>
</row>
</lyxtabular>
\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
<lyxtabular version="3" rows="5" columns="2">
<features tabularvalignment="middle" tabularwidth="100text%">
<column alignment="left" valignment="top" width="0pt">
<column alignment="right" valignment="top" width="8cm">
<row>
<cell alignment="left" valignment="top" usebox="none">
\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
</cell>
<cell alignment="right" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="left" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
External document ID:
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="left" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Document status:
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Draft
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="left" valignment="top" usebox="none">
\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
</cell>
<cell alignment="right" valignment="top" usebox="none">
\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
</cell>
</row>
<row>
<cell alignment="left" valignment="top" bottomline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" bottomline="true" usebox="none">
\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
</cell>
</row>
</lyxtabular>
\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<k<N_{\mathrm{DFT}}/2\\
X_{s}[k] & =\frac{\sqrt{2}}{N_{\mathrm{DFT}}}X[k] & \mathrm{for}\,\,\,k=N_{\mathrm{DFT}}/2
\end{eqnarray}
\end_inset
\end_layout
\begin_layout Subsection
Single sided power
\end_layout
\begin_layout Standard
The signal power can be computed from the single-sided amplitude spectrum
as
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\strikeout off
\uuline off
\uwave off
\noun off
\color none
\lang english
\begin_inset Formula $P(x[n])=\frac{1}{N_{\mathrm{DFT}}^{2}}\sum_{k=0}^{N_{\mathrm{DFT}}-1}X[k]X[k]^{*}$
\end_inset
\end_layout
\begin_layout Plain Layout
\lang english
-
\end_layout
\begin_layout Plain Layout
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\strikeout off
\uuline off
\uwave off
\noun off
\color none
\lang english
\begin_inset Formula $P(x[n])=\frac{2}{N_{\mathrm{DFT}}^{2}}\sum\limits _{k=0}^{N_{\mathrm{DFT}}/2}X[k]X[k]^{*}$
\end_inset
\end_layout
\begin_layout Plain Layout
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\strikeout off
\uuline off
\uwave off
\noun off
\color none
\lang english
\begin_inset Formula $P(x[n])=\frac{2}{N_{\mathrm{DFT}}^{2}}\sum\limits _{k=0}^{N_{\mathrm{DFT}}/2}\frac{N_{\mathrm{DFT}}}{2}X_{s}[k]\frac{N_{\mathrm{DFT}}}{2}X_{s}[k]^{*}=\frac{1}{2}\sum\limits _{k=0}^{N_{\mathrm{DFT}}/2}X_{s}[k]X_{s}[k]^{*}$
\end_inset
\end_layout
\end_inset
\begin_inset Formula
\begin{equation}
P(x[n])=\frac{1}{2}\sum_{k=0}^{N_{\mathrm{DFT}}/2}X_{s}[k]X_{s}[k]^{*},
\end{equation}
\end_inset
hence the signal power in frequency bin
\begin_inset Formula $k$
\end_inset
with with
\begin_inset Formula $\Delta f$
\end_inset
is
\begin_inset Formula
\begin{equation}
P_{k}=\frac{1}{2}\left\Vert X_{s}[k]\right\Vert ^{2}.
\end{equation}
\end_inset
The power spectral density for a single-sided spectrum is power per width
of the frequency bin:
\begin_inset Formula
\begin{equation}
\mathrm{PSD}_{k}=\frac{1}{\Delta f}P_{k}=\frac{N_{\mathrm{DFT}}}{f_{s}}P_{k}.
\end{equation}
\end_inset
\end_layout
\begin_layout Subsection
Welch' method of Cross power spectrum (CPS) estimation
\end_layout
\begin_layout Standard
Signals
\begin_inset Formula $x_{i}[n]$
\end_inset
.
Cross-power spectra as
\begin_inset Formula
\begin{equation}
C_{ij}=\frac{1}{2}X_{s,i}X_{s,j}^{*},
\end{equation}
\end_inset
where spectral averaging is used:
\begin_inset Formula
\begin{equation}
X_{s,i}X_{s,j}^{*}=\frac{1}{M}\sum_{m=0}^{M-1}\frac{\hat{X}_{s,i}^{(m)}\hat{X}_{s,j}^{(m)*}}{w_{p}},
\end{equation}
\end_inset
where
\begin_inset Formula $\hat{X}_{s,i}^{(m)}$
\end_inset
is the
\emph on
windowed
\emph default
time slot.
\end_layout
\begin_layout Section
Exponential time weighting on power spectra
\end_layout
\begin_layout Standard
Suppose we weigh exponential on
\begin_inset Formula $N$
\end_inset
power spectra, such that the total is
\begin_inset Formula
\begin{equation}
\boldsymbol{C}_{\mathrm{tot}}=\sum_{i=1}^{N}w_{i}\boldsymbol{C}_{i}
\end{equation}
\end_inset
Constraints:
\begin_inset Formula
\[
\sum_{i=1}^{N}w_{i}=1
\]
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
w_{i}=A\exp\left(-\alpha i\right)
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
w_{N}=0.01
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
\sum w_{i}=A\sum\exp\left(-\alpha i\right)
\end{equation}
\end_inset
\end_layout
\begin_layout Chapter
Impedance tube
\end_layout
\begin_layout Section
Design aspects
\end_layout
\begin_layout Itemize
Impedance tube should be constructed such that sound absorption measurements
can be done in compliance with the ASTM 1050 standard.
\end_layout
\begin_layout Itemize
The tube inner diameter is 50 mm, providing a cut-on frequency for higher
order modes of
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
\begin_inset Formula $f_{c}=\frac{c_{0}}{1.7D}=4047$
\end_inset
Hz
\end_layout
\begin_layout Plain Layout
According to standard:
\begin_inset Formula $f_{u}<Kc_{0}/D=0.586c_{0}/D=4031$
\end_inset
Hz.
\end_layout
\end_inset
4
\begin_inset space ~
\end_inset
kHz.
\end_layout
\begin_layout Itemize
The frequency range of interest for which the impedance tube will be used
is 40 Hz to 4 kHz, i.e.
two full octaves.
\end_layout
\begin_layout Itemize
A small close-able tube venting should be placed close to the speaker.
This venting should easily be opened and closed.
\end_layout
\begin_layout Itemize
A backing plate of at least 20 mm should be used to behind samples
\end_layout
\begin_layout Itemize
A rubber seal should be placed between the sound source and the tube to
avoid structure-borne sound:
\begin_inset Newline newline
\end_inset
\begin_inset Float figure
wide false
sideways false
status open
\begin_layout Plain Layout
\noindent
\align center
\begin_inset Graphics
filename img/radiator.png
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Caption Standard
\begin_layout Plain Layout
Fig.
copied from standard.
\end_layout
\end_inset
\end_layout
\end_inset
\end_layout
\begin_layout Subsection
Microphone distance
\end_layout
\begin_layout Itemize
Values computed for a sound speed of 344
\begin_inset space ~
\end_inset
m/s.
\end_layout
\begin_layout Itemize
The microphone distances calculated here are heart-heart distances.
Tube lengths should be calculated from these values based on the microphone
block geometry
\end_layout
\begin_layout Itemize
A length of 300 mm is used between source (speaker) and first microphone.
\end_layout
\begin_layout Itemize
For all frequencies of interest, we keep the microphone spacing
\begin_inset Formula $s$
\end_inset
such that:
\begin_inset Formula
\begin{equation}
0.1\pi\leq ks\leq0.8\pi
\end{equation}
\end_inset
Re-arranging: we can set the microphone spacing in terms of the upper and
lower frequencies:
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
\begin_inset Formula $0.1\pi\leq ks\leq0.8\pi$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $0.1\pi c_{0}\leq2\pi fs\leq0.8\pi c_{0}$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $0.05\frac{c_{0}}{f}\leq s\leq0.4\frac{c_{0}}{f}$
\end_inset
\end_layout
\end_inset
\begin_inset Formula
\begin{equation}
0.05\frac{c_{0}}{f_{l}}\leq s\leq0.4\frac{c_{0}}{f_{u}}
\end{equation}
\end_inset
We split up in three distances.
The first distance is used for the lowest frequency range, from 40 Hz and
upwards.
Setting
\begin_inset Formula $f_{l,1}$
\end_inset
to 40 Hz, and rounding of, we obtain
\begin_inset Formula $s_{1}=\boldsymbol{43}$
\end_inset
\begin_inset space ~
\end_inset
\series bold
cm
\series default
.
For distance
\begin_inset Formula $s_{1}$
\end_inset
, the upper measurement frequency is
\begin_inset Formula
\begin{equation}
f_{u,1}=0.4c_{0}/s_{1}=320\,\mathrm{Hz},
\end{equation}
\end_inset
which is exactly 8 times the lower frequency.
\begin_inset Newline newline
\end_inset
The smallest distance is chosen to set
\begin_inset Formula $f_{u,3}=4$
\end_inset
\begin_inset space ~
\end_inset
kHz.
Which gives a distance (rounded to millimeters!) of
\begin_inset Formula
\begin{equation}
s_{3}=0.4\frac{c_{0}}{f_{u,3}}=34.4\,\mathrm{mm}\Rightarrow\boldsymbol{34\,\mathrm{mm}}
\end{equation}
\end_inset
The lower measurement frequency for distance
\begin_inset Formula $s_{3}$
\end_inset
is:
\begin_inset Formula $f_{l,3}=506$
\end_inset
Hz.
Then, the distance
\begin_inset Formula $s_{2}$
\end_inset
is chosen to create an equal overlap on both the frequency range for distance
\begin_inset Formula $s_{1}$
\end_inset
as well as
\begin_inset Formula $s_{3}$
\end_inset
.
Hence, we set
\begin_inset Formula $f_{\mathrm{mid},2}=\left(f_{u,1}+f_{l,3}\right)/2=413$
\end_inset
\begin_inset space ~
\end_inset
Hz.
Hence, a bit rounded
\begin_inset Formula $s_{2}=\frac{0.4+0.05}{2}\frac{c_{0}}{f_{\mathrm{mid},2}}=19$
\end_inset
\begin_inset space ~
\end_inset
cm.
The following table summarizes the results:
\end_layout
\begin_layout Standard
\noindent
\align center
\begin_inset Tabular
<lyxtabular version="3" rows="4" columns="4">
<features booktabs="true" tabularvalignment="middle">
<column alignment="center" valignment="top">
<column alignment="right" valignment="top" width="0pt">
<column alignment="right" valignment="top">
<column alignment="right" valignment="top">
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Dimension
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Value
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Lower frequency
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
Upper frequency
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\begin_inset Formula $s_{1}$
\end_inset
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
430
\begin_inset space ~
\end_inset
mm
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
40
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
320
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\begin_inset Formula $s_{2}$
\end_inset
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
190
\begin_inset space ~
\end_inset
mm
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
90
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
724
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
\begin_inset Formula $s_{3}$
\end_inset
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
34
\begin_inset space ~
\end_inset
mm
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
491
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
<cell alignment="right" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text
\begin_layout Plain Layout
4047
\begin_inset space ~
\end_inset
Hz
\end_layout
\end_inset
</cell>
</row>
</lyxtabular>
\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}<n<N
\end{cases}
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
Hence:
\begin_inset Formula $\phi_{f}\equiv\phi_{N_{f}}$
\end_inset
is:
\begin_inset Formula
\begin{equation}
\phi_{/}=2\pi\Delta t\left[f_{L}N_{f}+\frac{1}{2}\left(N_{f}-1\right)\left(f_{U}-f_{L}\right)\right],
\end{equation}
\end_inset
and
\begin_inset Formula
\begin{equation}
\phi_{N}=\phi_{f}+2\pi\Delta t\left[\left(f_{U}+\varepsilon\right)N_{b}-\frac{1}{2}\left(N_{b}-1\right)\left(f_{U}+\varepsilon-f_{L}\right)\right],
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
and finally, after fulfilling the period, computing
\begin_inset Formula $K$
\end_inset
:
\begin_inset Formula
\begin{equation}
K=\left\lfloor \frac{\phi_{f}}{2\pi}+\Delta t\left[f_{U}N_{b}-\frac{1}{2}\left(N_{b}-1\right)\left(f_{U}-f_{L}\right)\right]\right\rfloor
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
Computing back the value of
\begin_inset Formula $\varepsilon$
\end_inset
:
\begin_inset Formula
\begin{equation}
\frac{\frac{1}{\Delta t}\left(K-\frac{\phi_{f}}{2\pi}\right)-f_{U}N_{b}+\frac{1}{2}\left(N_{b}-1\right)\left(f_{U}-f_{L}\right)}{\frac{1}{2}\left(N_{b}+1\right)}=\varepsilon.
\end{equation}
\end_inset
\end_layout
\begin_layout Subsection
Forward / backward exponential sweep
\end_layout
\begin_layout Standard
For a forward exponential sweep, the instantaneous frequency is defined
as:
\begin_inset Formula
\begin{equation}
f_{n}=f_{L}k^{\frac{n}{N}},
\end{equation}
\end_inset
where
\begin_inset Formula
\begin{equation}
k=\left(\frac{f_{U}+\varepsilon}{f_{L}}\right).
\end{equation}
\end_inset
Filling this in into the difference equation for the phase results in the
explicit solution for
\begin_inset Formula $\phi_{n}$
\end_inset
:
\begin_inset Formula
\begin{equation}
\phi_{n}=2\pi\Delta tf_{L}\frac{k^{\frac{n}{N}}-1}{\sqrt[N]{k}-1}.
\end{equation}
\end_inset
Tuning
\begin_inset Formula $K$
\end_inset
:
\begin_inset Formula
\begin{equation}
K=\left\lfloor \Delta tf_{L}\frac{k_{1}-1}{\sqrt[N]{k_{1}}-1}\right\rfloor ,
\end{equation}
\end_inset
The equation error
\begin_inset Formula $E$
\end_inset
is:
\begin_inset Formula
\begin{equation}
E=1-k+\frac{K}{\Delta tf_{L}}\left(\sqrt[N]{k}-1\right).
\end{equation}
\end_inset
\begin_inset Formula $E$
\end_inset
needs to be set to zero by adjusting
\begin_inset Formula $k$
\end_inset
.
This is a transcendental equation.
However, we already know that the starting value of
\begin_inset Formula $k$
\end_inset
is a quite proper starting value.
Therefore, we use Newton-Rhapson iterations:
\begin_inset Formula
\begin{equation}
\Delta k=\frac{-E}{\frac{\partial E}{\partial k}}=-\frac{E}{\frac{K}{\Delta tf_{L}}\frac{k^{\frac{1}{N}}}{Nk}}.
\end{equation}
\end_inset
\end_layout
\begin_layout Subsection
Continuous exponential sweep
\end_layout
\begin_layout Standard
As a starting value of
\begin_inset Formula $k$
\end_inset
, we set
\begin_inset Formula $k_{1}=\frac{f_{U}}{f_{L}}$
\end_inset
, then the phase after going forward in frequency is:
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
\phi_{f,1}=2\pi\Delta tf_{L}\frac{k_{1}-1}{\sqrt[N_{f}]{k}-1},
\end{equation}
\end_inset
where
\begin_inset Formula
\begin{equation}
k=\frac{f_{U}+\varepsilon}{f_{L}}\qquad;\qquad k_{1}=\frac{f_{U}}{f_{L}}
\end{equation}
\end_inset
such that
\begin_inset Formula
\begin{equation}
K=\left\lfloor \frac{\phi_{f,1}}{2\pi}+\Delta tf_{L}\frac{1-k_{1}}{\sqrt[N_{b}]{k_{1}^{-1}}-1}\right\rfloor ,
\end{equation}
\end_inset
The equation error for this case is:
\begin_inset Formula
\begin{equation}
E=\frac{k-1}{\sqrt[N_{f}]{k}-1}+k\frac{k^{-1}-1}{\sqrt[N_{b}]{1/k}-1}-\frac{K}{f_{L}\Delta t}.
\end{equation}
\end_inset
this equation is again solved numerically, using the starting value for
\begin_inset Formula $k$
\end_inset
as
\begin_inset Formula $k_{1}$
\end_inset
.
The derivative of
\begin_inset Formula $E$
\end_inset
to
\begin_inset Formula $k$
\end_inset
is:
\begin_inset Formula
\begin{equation}
\frac{\mathrm{d}E}{\mathrm{d}k}=\frac{\frac{1}{k}-1}{k^{-\frac{1}{N_{b}}}-1}+\frac{1}{k^{\frac{1}{N_{f}}}-1}+\frac{1}{k\left(k^{-\frac{1}{N_{b}}}-1\right)}+\frac{k^{\frac{1}{N_{f}}}\left(1-k\right)}{N_{f}k\left(k^{\frac{1}{N_{f}}}-1\right)^{2}}+\frac{k^{-\frac{1}{N_{b}}}\left(-1+\frac{1}{k}\right)}{N_{b}\left(k^{-\frac{1}{N_{b}}}-1\right)^{2}}
\end{equation}
\end_inset
\end_layout
\begin_layout Section
Transfer function estimation using cross-power spectra
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
C_{ij}=\frac{1}{2}P_{i}P_{j}^{*}
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
Estimation of
\begin_inset Formula
\begin{equation}
G_{ij}=\frac{P_{j}}{P_{i}}\approx\frac{C_{ji}}{C_{ii}}
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
Or:
\begin_inset Formula
\begin{equation}
G_{ij}=\frac{P_{j}}{P_{i}}=\frac{C_{jj}}{C_{ij}}
\end{equation}
\end_inset
\end_layout
\begin_layout Section
Impulse response and transfer function estimation using coherent periodic
averaging
\end_layout
\begin_layout Standard
Input signal:
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
x[n]=x[n+N]\qquad N\in\mathbb{Z}
\end{equation}
\end_inset
Generating a periodic sequence, avoid temporal aliasing.
\end_layout
\begin_layout Standard
We would like to know the transfer function
\begin_inset Formula $h[t]$
\end_inset
, which is defined such that the output
\begin_inset Formula $y[n]$
\end_inset
is:
\begin_inset Formula
\begin{equation}
y[n]=h*x[n]=\sum_{m=-\infty}^{\infty}h[n-m]x[n]
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
Suppose the input signal is a periodic sine sweep, with a quiescent tail.
We can define this sweep as the convolved
\begin_inset Formula
\begin{equation}
x[n]=s[n]*\delta[n-m]
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
The output is then:
\begin_inset Formula
\begin{equation}
y[n]=h*s*\delta+e,
\end{equation}
\end_inset
where
\begin_inset Formula $e$
\end_inset
is the uncorrelated noise.
Then by pre-multiplying with
\begin_inset Formula $s^{-1}$
\end_inset
, i.e.
the inverse of
\begin_inset Formula $s$
\end_inset
, this results in:
\begin_inset Formula
\begin{equation}
s^{-1}*y[n]=s^{-1}*h*s*\delta[n-m]+s^{-1}*e[n].
\end{equation}
\end_inset
As the convolution is associative
\begin_inset Formula $s^{-1}*s$
\end_inset
cancel out, and this can be written as:
\begin_inset Formula
\begin{equation}
s^{-1}*y[n]=h*\delta[n-m]+s^{-1}*e[n],
\end{equation}
\end_inset
hence the impulse response can be written as
\begin_inset Formula
\begin{equation}
h=s^{-1}*\left(y[n]-e[n]\right)
\end{equation}
\end_inset
\end_layout
\begin_layout Section
Transfer function estimation using least squares estimation and a FIR filter
\end_layout
\begin_layout Standard
This is a time domain method to estimate the transfer function by minimizing
the error.
For a FIR system we can write
\begin_inset Formula
\begin{equation}
y[n]=\sum_{l=0}^{L-1}w[n-l]x[n-l]+d[n]=\boldsymbol{w}\cdot\boldsymbol{x}[n]+d[n],
\end{equation}
\end_inset
where
\begin_inset Formula $\boldsymbol{w}$
\end_inset
is the vector of FIR coefficients
\begin_inset Formula $d[n]$
\end_inset
is an unknown disturbance (noise), and
\begin_inset Formula $\boldsymbol{x}[n]$
\end_inset
the time history at discrete time point
\begin_inset Formula $n$
\end_inset
containing
\begin_inset Formula $L$
\end_inset
points.
If we write
\begin_inset Formula
\begin{equation}
e[n]=y[n]-\boldsymbol{w}\cdot\boldsymbol{x}[n],
\end{equation}
\end_inset
we can generate the cost function
\begin_inset Formula
\begin{equation}
J=\sum_{n=0}^{N-1}e^{2}[n]
\end{equation}
\end_inset
to be minimized, we can estimate the optimal filter
\begin_inset Formula $\boldsymbol{w}$
\end_inset
by taking the derivative of
\begin_inset Formula $J$
\end_inset
w.r.t.
to
\begin_inset Formula $\boldsymbol{w}$
\end_inset
and setting that to 0.
This results in:
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
\begin_inset Formula $\frac{\partial J}{\partial\boldsymbol{w}}=\frac{\partial}{\partial\boldsymbol{w}}\sum_{n=0}^{N-1}\left[e[n]^{2}\right]=\sum_{n=0}^{N-1}\left[2e[n]\frac{\partial e[n]}{\partial\boldsymbol{w}}\right]$
\end_inset
\end_layout
\begin_layout Plain Layout
Working out:
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $\frac{\partial J}{\partial\boldsymbol{w}}=\sum_{n=0}^{N-1}\left[2e[n]\boldsymbol{x}[n]\right]$
\end_inset
\end_layout
\begin_layout Plain Layout
Setting this to
\begin_inset Formula $\boldsymbol{0}$
\end_inset
:
\end_layout
\begin_layout Plain Layout
\begin_inset Formula
\[
\sum_{n=0}^{N-1}\left[2\left(y[n]-\boldsymbol{w}\cdot\boldsymbol{x}[n]\right)\boldsymbol{x}[n]\right]=\boldsymbol{0}
\]
\end_inset
\end_layout
\begin_layout Plain Layout
Defining the matrix
\begin_inset Formula $X_{nl}$
\end_inset
containing the time history at time at time instance
\begin_inset Formula $n$
\end_inset
and history point
\begin_inset Formula $l$
\end_inset
, and
\begin_inset Formula $Y_{nl}$
\end_inset
containing the output signal at time instance
\begin_inset Formula $n$
\end_inset
, we can write
\begin_inset Formula $\hat{\boldsymbol{y}}=\boldsymbol{X}\cdot\boldsymbol{w}$
\end_inset
, such that:
\begin_inset Formula
\[
\sum_{n=0}^{N-1}\left[2\left(y[n]-\boldsymbol{w}\cdot\boldsymbol{x}[n]\right)\boldsymbol{x}[n]\right]=\boldsymbol{0}
\]
\end_inset
\end_layout
\begin_layout Plain Layout
For all
\end_layout
\end_inset
\end_layout
\begin_layout Section
Vector fitting
\end_layout
\begin_layout Itemize
Method to estimate poles and zeros (rational transfer function estimation)
of an arbitrary transfer function.
Guaranteed stable poles.
\end_layout
\begin_layout Standard
Suppose we want to fit the transfer function
\begin_inset Formula $f(s)$
\end_inset
using a rational function approximation as:
\begin_inset Formula
\begin{equation}
f(s)\approx\sum_{n=1}^{N}\frac{c_{n}}{s-a_{n}}+d+sh,
\end{equation}
\end_inset
where
\begin_inset Formula $a_{n}$
\end_inset
are the poles.
We introduce an
\emph on
unknown
\emph default
function
\begin_inset Formula $\sigma(s)$
\end_inset
which shares the same poles as
\begin_inset Formula $f(s)$
\end_inset
, but with different zeros:
\begin_inset Formula
\begin{equation}
\sigma(s)=\sum_{n=1}^{N}\frac{\tilde{c}_{n}}{s-a_{n}}+1
\end{equation}
\end_inset
We create the augmented problem:
\begin_inset Formula
\begin{equation}
\left[\begin{array}{c}
\sigma(s)f(s)\\
\sigma(s)
\end{array}\right]\approx\left[\begin{array}{c}
\sum\limits _{n=1}^{N}\frac{c_{n}}{s-a_{n}}+d+sh\\
\sum_{n=1}^{N}\frac{\tilde{c}_{n}}{s-a_{n}}+k(s)
\end{array}\right].
\end{equation}
\end_inset
So we assume that we can fit
\begin_inset Formula $\sigma(s)f(s)$
\end_inset
with the top equation, and
\begin_inset Formula $\sigma(s)$
\end_inset
with the bottom equation.
By multiplying the bottom equation with
\begin_inset Formula $f(s)$
\end_inset
, we find:
\begin_inset Formula
\begin{equation}
\underbrace{\sum\limits _{n=1}^{N}\frac{c_{n}}{s-a_{n}}+d+sh}_{\left(\sigma f\right)_{\mathrm{fit}}}=\underbrace{\left(\sum_{n=1}^{N}\frac{\tilde{c}_{n}}{s-a_{n}}+1\right)f(s)}_{\sigma_{\mathrm{fit}}(s)f(s)},\label{eq:sigmaf_eq_sigma_f}
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
For a certain set of starting poles
\begin_inset Formula $\overline{a}_{n}$
\end_inset
, we are able to fit the zeros.
Eq.
\begin_inset space ~
\end_inset
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:sigmaf_eq_sigma_f"
\end_inset
is a linear least squares problem, which is used to fit
\end_layout
\begin_layout Section
Soft gain transitioning
\end_layout
\begin_layout Standard
Equation governing gain evolution, running at each sample:
\begin_inset Formula
\begin{equation}
\Delta g=\alpha\left(g_{\mathrm{required}}-g_{\mathrm{old}}\right),
\end{equation}
\end_inset
which can be written as:
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
\begin_inset Formula $G[z]\left(1-z^{-1}\right)=\alpha\left(G_{\mathrm{required}}[z]-z^{-1}G[z]\right)$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $G[z]\left(1-z^{-1}\right)+\alpha z^{-1}G[z]=\alpha G_{\mathrm{required}}[z]$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $G[z]=\frac{\alpha}{1+\left(\alpha-1\right)z^{-1}}G_{\mathrm{required}}[z]$
\end_inset
\end_layout
\end_inset
\begin_inset Formula
\begin{equation}
G[z]=\frac{\alpha}{1+\left(\alpha-1\right)z^{-1}}G_{\mathrm{required}}(z),
\end{equation}
\end_inset
which is an approximate first order digital low-pass filter:
\begin_inset Note Note
status collapsed
\begin_layout Plain Layout
Filling in for
\begin_inset Formula $H(s)=H[z=e^{sT}],$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $G\left(s\right)=\frac{\alpha}{1+\left(\alpha-1\right)\exp\left(-s/f_{s}\right)}G_{\mathrm{required}}(s)$
\end_inset
,for
\begin_inset Formula $s\ll f_{s}$
\end_inset
can be written as:
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $G\left(s\right)\approx\frac{\alpha}{\alpha+\left(1-\alpha\right)s/f_{s}}G_{\mathrm{required}}(s)$
\end_inset
which has a pole at approximately:
\begin_inset Formula $s_{p}=f_{s}\left[1-\frac{1}{1-\alpha}\right]$
\end_inset
\end_layout
\begin_layout Plain Layout
Fill in:
\begin_inset Formula $\alpha=0.1\Rightarrow s_{p}=f_{s}\left[1-\frac{1}{0.9}\right]$
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Formula $s=f_{s}\left[-9\right]$
\end_inset
\end_layout
\end_inset
\begin_inset Formula
\begin{equation}
G\left(s\right)\approx\frac{\alpha}{\alpha+s/f_{s}}G_{\mathrm{required}}(s)\qquad s\ll f_{s},0<\alpha\ll1
\end{equation}
\end_inset
Hence for the transition time constant,
\begin_inset Formula $s_{p}=-\frac{2\pi}{\tau_{p}}$
\end_inset
, we find:
\begin_inset Note Note
status open
\begin_layout Plain Layout
\begin_inset Formula $\alpha=1-\frac{1}{1+\frac{2\pi}{\tau_{p}f_{s}}}$
\end_inset
\end_layout
\end_inset
\begin_inset Formula
\begin{equation}
\alpha=1-\frac{1}{1+\frac{2\pi}{\tau_{p}f_{s}}}
\end{equation}
\end_inset
\end_layout
\begin_layout Chapter
Programming aspects
\end_layout
\begin_layout Section
Array ordering
\end_layout
\begin_layout Itemize
A 2D array of is called in
\emph on
row-major order
\emph default
if each of the rows is subsequentially put in memory.
This means, we can find element
\begin_inset Formula $A_{\mathrm{row},\mathrm{col}}=A_{ij}$
\end_inset
in memory
\begin_inset Formula $A[i+j\mathrm{n_{rows}}]$
\end_inset
.
This is also called a C-contiguous array.
\end_layout
\begin_layout Itemize
A 2D array of is called in
\emph on
column-major order
\emph default
if each of the columns is subsequentially put in memory.
This means, we can find element
\begin_inset Formula $A_{\mathrm{row},\mathrm{col}}=A_{ij}$
\end_inset
in memory
\begin_inset Formula $A[in_{\mathrm{cols}}+j]$
\end_inset
.
This is also called a Fortran-contiguous array.
For the
\family typewriter
dmat
\family default
structure in
\family typewriter
LASP
\family default
, we have chosen to use the Fortran-contiguous convention.
\end_layout
\begin_layout Standard
Numpy works under the hood with
\family typewriter
strides
\family default
, strides are the values with which each array index needs to be multiplied
to obtain the value of each index in an
\begin_inset Formula $N$
\end_inset
-dimensional array.
\family typewriter
Strides
\family default
is a tuple containing
\begin_inset Formula $(a,b,c,\dots$
\end_inset
).
Such that
\begin_inset Formula $A[i,j,k,\dots]=A[ai+bj+ck+\dots]$
\end_inset
.
Strides are a method to have an arbitrary memory layout, and also provide
a very easy way to transpose an array over one dimension (just swap the
stride values in the tuple).
\end_layout
\begin_layout Chapter
Beamforming
\end_layout
\begin_layout Section
Near-field cardioid beamformer
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
\boldsymbol{x}(\omega)=\boldsymbol{a}(\omega)s,
\end{equation}
\end_inset
where
\begin_inset Formula $s$
\end_inset
is the source signal, and
\begin_inset Formula $\boldsymbol{a}(\omega)$
\end_inset
the transfer from the source signal to the measured signal
\begin_inset Formula $\boldsymbol{x}$
\end_inset
.
Now, we would like to generate a signal
\begin_inset Formula $y=\boldsymbol{w}^{H}\boldsymbol{x}$
\end_inset
, that is sensitive to signals at position
\begin_inset Formula $\boldsymbol{r}_{s}$
\end_inset
and not sensitive to signals at other positions.
Therefore, the optimal beamformer is such that:
\begin_inset Formula
\begin{align}
\boldsymbol{w}^{H}\boldsymbol{a}s_{\boldsymbol{r}=\boldsymbol{r}_{0}} & =1,\\
\boldsymbol{w}^{H}\boldsymbol{a}s_{\boldsymbol{r}\neq\boldsymbol{r}_{0}} & =0.
\end{align}
\end_inset
Now, the sensitivity to a source is:
\begin_inset Formula
\begin{equation}
\boldsymbol{a}=\left\{ \begin{array}{c}
\exp\left(-i\omega R_{1}/c_{0}\right)/R_{1}\\
\exp\left(-i\omega R_{2}/c_{0}\right)/R_{2}\\
\vdots\\
\\
\end{array}\right\}
\end{equation}
\end_inset
\end_layout
\begin_layout Standard
\begin_inset CommandInset bibtex
LatexCommand bibtex
btprint "btPrintCited"
bibfiles "lasp"
options "bibtotoc,plain"
\end_inset
\end_layout
\end_body
\end_document
Reference in New Issue View Git Blame Copy Permalink