1199 lines
45 KiB
HTML
1199 lines
45 KiB
HTML
<?xml version="1.0" encoding="utf-8"?>
|
||
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
|
||
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
|
||
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
|
||
<head>
|
||
<!-- 2021-04-26 lun. 17:04 -->
|
||
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
|
||
<title>Spectral Analysis using Matlab</title>
|
||
<meta name="author" content="Dehaeze Thomas" />
|
||
<meta name="generator" content="Org Mode" />
|
||
<link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
|
||
<script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
|
||
<script>
|
||
MathJax = {
|
||
svg: {
|
||
scale: 1,
|
||
fontCache: "global"
|
||
},
|
||
tex: {
|
||
tags: "ams",
|
||
multlineWidth: "%MULTLINEWIDTH",
|
||
tagSide: "right",
|
||
macros: {bm: ["\\boldsymbol{#1}",1],},
|
||
tagIndent: ".8em"
|
||
}
|
||
};
|
||
</script>
|
||
<script id="MathJax-script" async
|
||
src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js"></script>
|
||
</head>
|
||
<body>
|
||
<div id="org-div-home-and-up">
|
||
<a accesskey="h" href="../index.html"> UP </a>
|
||
|
|
||
<a accesskey="H" href="../index.html"> HOME </a>
|
||
</div><div id="content">
|
||
<h1 class="title">Spectral Analysis using Matlab</h1>
|
||
<div id="table-of-contents">
|
||
<h2>Table of Contents</h2>
|
||
<div id="text-table-of-contents">
|
||
<ul>
|
||
<li><a href="#org53dcf97">1. Spectral Analysis - Basics</a>
|
||
<ul>
|
||
<li><a href="#org3ed291d">1.1. Sensitivity of the instrumentation</a></li>
|
||
<li><a href="#org91e6c0f">1.2. Convert the time domain from volts to velocity</a></li>
|
||
<li><a href="#org80959fd">1.3. Power Spectral Density and Amplitude Spectral Density</a></li>
|
||
<li><a href="#org07905ac">1.4. Modification of a signal’s Power Spectral Density when going through an LTI system</a></li>
|
||
<li><a href="#orgdb7af23">1.5. From PSD of the velocity to the PSD of the displacement</a></li>
|
||
<li><a href="#orga382ecf">1.6. Cumulative Power/Amplitude Spectrum</a></li>
|
||
</ul>
|
||
</li>
|
||
<li><a href="#org7953d16">2. Time domain signal that approximate a PSD - TF technique</a>
|
||
<ul>
|
||
<li><a href="#org4a1b3f8">2.1. Signal’s PSD</a></li>
|
||
<li><a href="#orgeba7bdf">2.2. Transfer Function that approximate the ASD</a></li>
|
||
<li><a href="#orgbcffeed">2.3. Generated Time domain signal</a></li>
|
||
<li><a href="#org676e750">2.4. Comparison of the Power Spectral Densities</a></li>
|
||
<li><a href="#org8669bb0">2.5. Simulink</a></li>
|
||
</ul>
|
||
</li>
|
||
<li><a href="#org7073f62">3. Time domain signal that approximate a PSD - IFFT technique</a>
|
||
<ul>
|
||
<li><a href="#orge49d520">3.1. Signal’s PSD</a></li>
|
||
<li><a href="#orgcb7dbc6">3.2. Algorithm</a></li>
|
||
<li><a href="#orgf14004b">3.3. Obtained Time Domain Signal</a></li>
|
||
<li><a href="#org5d770d3">3.4. PSD Comparison</a></li>
|
||
</ul>
|
||
</li>
|
||
<li><a href="#org9515c46">4. Compute the Noise level and Signal level from PSD</a>
|
||
<ul>
|
||
<li><a href="#org073b90d">4.1. Time Domain Signal</a></li>
|
||
<li><a href="#org86ba021">4.2. Estimation of the magnitude of a deterministic signal</a></li>
|
||
<li><a href="#org833e18d">4.3. Estimation of the noise level</a></li>
|
||
</ul>
|
||
</li>
|
||
<li><a href="#orgd325fcb">5. Further Notes</a>
|
||
<ul>
|
||
<li><a href="#orgf74d5dc">5.1. PSD of ADC quantization noise</a></li>
|
||
</ul>
|
||
</li>
|
||
</ul>
|
||
</div>
|
||
</div>
|
||
|
||
<p>
|
||
This document presents the mathematics as well as the matlab scripts to do various spectral analysis on a measured signal.
|
||
</p>
|
||
|
||
<p>
|
||
Some matlab documentation about Spectral Analysis can be found <a href="https://fr.mathworks.com/help/signal/ug/spectral-analysis.html">here</a>.
|
||
</p>
|
||
|
||
<p>
|
||
First, in section <a href="#org9f5a786">1</a>, some basics of spectral analysis are presented.
|
||
</p>
|
||
|
||
<p>
|
||
In some cases, we want to generate a time domain signal with defined Power Spectral Density.
|
||
Two methods are presented in sections <a href="#org03ec73d">2</a> and <a href="#org401f4ee">3</a>.
|
||
</p>
|
||
|
||
<p>
|
||
Finally, some notes are done on how to compute the noise level and signal level from a given Power Spectral Density in section <a href="#orga227890">4</a>.
|
||
</p>
|
||
|
||
<div id="outline-container-org53dcf97" class="outline-2">
|
||
<h2 id="org53dcf97"><span class="section-number-2">1</span> Spectral Analysis - Basics</h2>
|
||
<div class="outline-text-2" id="text-1">
|
||
<p>
|
||
<a id="org9f5a786"></a>
|
||
</p>
|
||
<p>
|
||
In this section, the basics of spectral analysis is presented with the associated Matlab commands.
|
||
</p>
|
||
|
||
<p>
|
||
This include:
|
||
</p>
|
||
<ul class="org-ul">
|
||
<li>how to compute the Power Spectral Density (PSD) and the Amplitude Spectral Density (ASD)</li>
|
||
<li>how to take into account the sensitivity of the sensor and of the electronics</li>
|
||
<li>how to compute the Cumulative Power Spectrum (CPS) and the Cumulative Amplitude Spectrum (CAS)</li>
|
||
</ul>
|
||
</div>
|
||
|
||
<div id="outline-container-org3ed291d" class="outline-3">
|
||
<h3 id="org3ed291d"><span class="section-number-3">1.1</span> Sensitivity of the instrumentation</h3>
|
||
<div class="outline-text-3" id="text-1-1">
|
||
<p>
|
||
A typical measurement setup is shown in figure <a href="#org7b81fb1">1</a> where we measure a physical signal which is here a velocity \(v(t)\) using a geophone.
|
||
The geophone has some dynamics that we represent with \(G_g(s)\), its output a voltage.
|
||
The output of the geophone is then amplified by a voltage amplifier with a transfer function \(G_a(s)\).
|
||
</p>
|
||
|
||
<p>
|
||
Finally, the signal is discretised with an ADC and we obtain \(x(k\Delta t)\).
|
||
</p>
|
||
|
||
<p>
|
||
To obtain the real physical quantity \(v(t)\) as measured by the sensor from the obtained signal \(x(k\Delta t)\), one have to know the sensitivity of the sensors and electronics used (\(G_g(s)\) and \(G_a(s)\)).
|
||
</p>
|
||
|
||
|
||
<div id="org7b81fb1" class="figure">
|
||
<p><img src="figs/velocity_to_voltage.png" alt="velocity_to_voltage.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 1: </span>Schematic of the instrumentation used for the measurement</p>
|
||
</div>
|
||
|
||
<p>
|
||
Let’s say, we know that the sensitivity of the geophone used is
|
||
\[ G_g(s) = G_0 \frac{\frac{s}{2\pi f_0}}{1 + \frac{s}{2\pi f_0}} \quad \left[\frac{V}{m/s}\right] \]
|
||
</p>
|
||
|
||
<p>
|
||
The parameters are defined below.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> G0 = 88; <span class="org-comment">% Sensitivity [V/(m/s)]</span>
|
||
f0 = 2; <span class="org-comment">% Cut-off frequency [Hz]</span>
|
||
|
||
Gg = G0<span class="org-type">*</span>(s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>f0)<span class="org-type">/</span>(1<span class="org-type">+</span>s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>f0);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
And the dynamics of the amplifier in the bandwidth of interest is just a gain: \(G_a(s) = 1000\).
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> Ga = 1000;
|
||
</pre>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org91e6c0f" class="outline-3">
|
||
<h3 id="org91e6c0f"><span class="section-number-3">1.2</span> Convert the time domain from volts to velocity</h3>
|
||
<div class="outline-text-3" id="text-1-2">
|
||
<p>
|
||
Let’s here try to obtain the time domain signal \(v(t)\) from the measurement \(x\).
|
||
</p>
|
||
|
||
<p>
|
||
If \({G_a(s)}^{-1} {G_g(s)}^{-1}\) is proper, we can simulate this dynamical system to go from the voltage \(x\) to the velocity \(v\) as shown in figure <a href="#orgd072aa0">2</a>.
|
||
</p>
|
||
|
||
<p>
|
||
If \({G_a(s)}^{-1} {G_g(s)}^{-1}\) is not proper, we add low pass filters at high frequency to make the system proper.
|
||
</p>
|
||
|
||
|
||
<div id="orgd072aa0" class="figure">
|
||
<p><img src="figs/voltage_to_velocity.png" alt="voltage_to_velocity.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 2: </span>Schematic of the instrumentation used for the measurement</p>
|
||
</div>
|
||
|
||
|
||
<p>
|
||
Let’s load the measured \(x\).
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> data = load(<span class="org-string">'mat/data_028.mat'</span>, <span class="org-string">'data'</span>); data = data.data;
|
||
|
||
t = data(<span class="org-type">:</span>, 3); <span class="org-comment">% Time vector [s]</span>
|
||
x = data(<span class="org-type">:</span>, 1)<span class="org-type">-</span>mean(data(<span class="org-type">:</span>, 1)); <span class="org-comment">% The offset if removed (coming from the voltage amplifier) [v]</span>
|
||
|
||
dt = t(2)<span class="org-type">-</span>t(1); <span class="org-comment">% Sampling Time [s]</span>
|
||
Fs = 1<span class="org-type">/</span>dt; <span class="org-comment">% Sampling Frequency [Hz]</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We simulate this system with matlab using the <code>lsim</code> command.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> v = lsim(inv(Gg<span class="org-type">*</span>Ga), x, t);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
And we plot the obtained velocity
|
||
</p>
|
||
|
||
<div id="orga9c51ca" class="figure">
|
||
<p><img src="figs/velocity_time.png" alt="velocity_time.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 3: </span>Computed Velocity from the measured Voltage</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org80959fd" class="outline-3">
|
||
<h3 id="org80959fd"><span class="section-number-3">1.3</span> Power Spectral Density and Amplitude Spectral Density</h3>
|
||
<div class="outline-text-3" id="text-1-3">
|
||
<p>
|
||
From the Matlab documentation:
|
||
</p>
|
||
<blockquote>
|
||
<p>
|
||
The goal of spectral estimation is to describe the distribution (over frequency) of the power contained in a signal, based on a finite set of data.
|
||
</p>
|
||
</blockquote>
|
||
|
||
<p>
|
||
We now have the velocity \(v(t)\ [m/s]\) in the time domain.
|
||
</p>
|
||
|
||
<p>
|
||
The Power Spectral Density (PSD) \(S_v(f)\) of the time domain \(v(t)\) can be computed using the following equation:
|
||
\[ S_v(f) = \frac{1}{f_s} \sum_{m=-\infty}^{\infty} R_{xx}(m) e^{-j 2 \pi m f / f_s} \ \left[\frac{(m/s)^2}{Hz}\right] \]
|
||
where
|
||
</p>
|
||
<ul class="org-ul">
|
||
<li>\(f_s\) is the sampling frequency in Hz</li>
|
||
<li>\(R_{xx}\) is the autocorrelation</li>
|
||
</ul>
|
||
|
||
<p>
|
||
The PSD represents the distribution of the (average) signal power over frequency.
|
||
\(S_v(f)df\) is the infinitesimal power in the band \((f-df/2, f+df/2)\) and the total power in the signal is obtained by integrating these infinitesimal contributions.
|
||
</p>
|
||
|
||
|
||
<p>
|
||
To compute the Power Spectral Density with matlab, we use the <code>pwelch</code> function (<a href="https://fr.mathworks.com/help/signal/ref/pwelch.html?s_tid=doc_ta">documentation</a>).
|
||
The use of the <code>pwelch</code> function is:
|
||
<code>[pxx,w] = pwelch(x,window,noverlap,nfft, fs)</code>
|
||
with:
|
||
</p>
|
||
<ul class="org-ul">
|
||
<li><code>x</code> is the discrete time signal</li>
|
||
<li><code>window</code> is a window that is used to smooth the obtained PSD</li>
|
||
<li><code>overlap</code> can be used to have some overlap from section to section</li>
|
||
<li><code>nfft</code> specifies the number of FFT points for the PSD</li>
|
||
<li><code>fs</code> is the sampling frequency of the data <code>x</code> in Hertz</li>
|
||
</ul>
|
||
|
||
<p>
|
||
As explained in (<a href="#citeproc_bib_item_2">Schmid 2012</a>), it is recommended to use the <code>pwelch</code> function the following way:
|
||
</p>
|
||
|
||
<p>
|
||
First, define a window (preferably the <code>hanning</code> one) by specifying the averaging factor <code>na</code>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(v);
|
||
na = 8;
|
||
win = hanning(floor(nx<span class="org-type">/</span>na));
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Then, compute the power spectral density \(S_v\) and the associated frequency vector \(f\) with no overlap.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> [Sv, f] = pwelch(v, win, 0, [], Fs);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The obtained PSD is shown in figure <a href="#org55af140">4</a>.
|
||
</p>
|
||
|
||
|
||
<div id="org55af140" class="figure">
|
||
<p><img src="figs/psd_velocity.png" alt="psd_velocity.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 4: </span>Power Spectral Density of the measured velocity</p>
|
||
</div>
|
||
|
||
<p>
|
||
The Amplitude Spectral Density (ASD) is defined as the square root of the Power Spectral Density and is shown in figure <a href="#org1e82a61">5</a>.
|
||
</p>
|
||
\begin{equation}
|
||
\Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right]
|
||
\end{equation}
|
||
|
||
|
||
<div id="org1e82a61" class="figure">
|
||
<p><img src="figs/asd_velocity.png" alt="asd_velocity.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 5: </span>Power Spectral Density of the measured velocity</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org07905ac" class="outline-3">
|
||
<h3 id="org07905ac"><span class="section-number-3">1.4</span> Modification of a signal’s Power Spectral Density when going through an LTI system</h3>
|
||
<div class="outline-text-3" id="text-1-4">
|
||
<p>
|
||
Instead of computing the time domain velocity before computing the Power Spectral Density, we could have directly computed the PSD of the measured voltage \(x\) and then take into account the sensitivity of the measurement devices to have the PSD of the velocity.
|
||
</p>
|
||
|
||
<p>
|
||
To do so, we use the fact that a signal \(u\) with a PSD \(S_{uu}\) going through a LTI system \(G_(s)\) (figure <a href="#org5d3e0cf">6</a>) will generate a signal \(y\) with a PSD:
|
||
</p>
|
||
\begin{equation}
|
||
S_{yy}(\omega) = \left|G(j\omega)\right|^2 S_{uu}(\omega)
|
||
\end{equation}
|
||
|
||
|
||
<div id="org5d3e0cf" class="figure">
|
||
<p><img src="figs/velocity_to_voltage_psd.png" alt="velocity_to_voltage_psd.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 6: </span>Schematic of the instrumentation used for the measurement</p>
|
||
</div>
|
||
|
||
<p>
|
||
Similarly, the ASD of \(y\) is:
|
||
</p>
|
||
\begin{equation}
|
||
\Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{uu}(\omega)
|
||
\end{equation}
|
||
|
||
<p>
|
||
Thus, we could have computed the PSD of \(x\) and then obtain the PSD of the velocity with:
|
||
\[ S_{v}(\omega) = |G_a(j\omega) G_g(j\omega)|^{-1} S_{x}(\omega) \]
|
||
</p>
|
||
|
||
<p>
|
||
The PSD of \(x\) is computed below.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(x);
|
||
na = 8;
|
||
win = hanning(floor(nx<span class="org-type">/</span>na));
|
||
[Sx, f] = pwelch(x, win, 0, [], Fs);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
And the PSD of \(v\) is obtained with the below code.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> Svv = Sx<span class="org-type">.*</span>squeeze(abs(freqresp(inv(Gg<span class="org-type">*</span>Ga), f, <span class="org-string">'Hz'</span>)))<span class="org-type">.^</span>2;
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The result is compare with the PSD computed from the \(v\) signal obtained with the <code>lsim</code> command in figure <a href="#org39f0907">7</a>.
|
||
</p>
|
||
|
||
|
||
<div id="org39f0907" class="figure">
|
||
<p><img src="figs/psd_velocity_lti_method.png" alt="psd_velocity_lti_method.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 7: </span>Obtain PSD of the velocity using the formula (<a href="./figs/psd_velocity_lti_method.png">png</a>, <a href="./figs/psd_velocity_lti_method.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgdb7af23" class="outline-3">
|
||
<h3 id="orgdb7af23"><span class="section-number-3">1.5</span> From PSD of the velocity to the PSD of the displacement</h3>
|
||
<div class="outline-text-3" id="text-1-5">
|
||
<p>
|
||
Similarly to what has been done in the last section, we can consider the displacement \(d\) can be obtained from the velocity \(v\) by going through an LTI system \(1/s\) as shown in figure <a href="#orgeae869f">8</a>.
|
||
</p>
|
||
|
||
|
||
<div id="orgeae869f" class="figure">
|
||
<p><img src="figs/velocity_to_displacement_psd.png" alt="velocity_to_displacement_psd.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 8: </span>Schematic of the instrumentation used for the measurement</p>
|
||
</div>
|
||
|
||
<p>
|
||
We then have the relation between the PSD of \(d\) and the PSD of \(v\):
|
||
</p>
|
||
\begin{equation}
|
||
S_{dd}(\omega) = \left|\frac{1}{j \omega}\right|^2 S_{vv}(\omega)
|
||
\end{equation}
|
||
|
||
<p>
|
||
Using a frequency variable \(f\) in Hz:
|
||
</p>
|
||
\begin{equation}
|
||
S_{dd}(f) = \left| \frac{1}{j 2\pi f} \right|^2 S_{vv}(f)
|
||
\end{equation}
|
||
|
||
<p>
|
||
For the Amplitude Spectral Density:
|
||
</p>
|
||
\begin{equation}
|
||
\Gamma_{dd}(f) = \frac{1}{2\pi f} \Gamma_{vv}(f)
|
||
\end{equation}
|
||
|
||
<p>
|
||
Note here that the PSD (and ASD) of one variable and its derivatives/integrals are equal at one particular frequency \(f = 1\ rad/s \approx 0.16\ Hz\):
|
||
</p>
|
||
\begin{equation}
|
||
S_{xx}(\omega = 1) = S_{vv}(\omega = 1)
|
||
\end{equation}
|
||
|
||
|
||
<p>
|
||
With Matlab, the PSD of the displacement can be computed from the PSD of the velocity with the following code.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> Sd = Sv<span class="org-type">.*</span>(1<span class="org-type">./</span>(2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>f))<span class="org-type">.^</span>2;
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The obtained PSD of the displacement can be seen in figure <a href="#orgfd4a6af">9</a>.
|
||
</p>
|
||
|
||
<div id="orgfd4a6af" class="figure">
|
||
<p><img src="figs/psd_velocity_displacement.png" alt="psd_velocity_displacement.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 9: </span>PSD of the Velocity and Displacement (<a href="./figs/psd_velocity_displacement.png">png</a>, <a href="./figs/psd_velocity_displacement.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orga382ecf" class="outline-3">
|
||
<h3 id="orga382ecf"><span class="section-number-3">1.6</span> Cumulative Power/Amplitude Spectrum</h3>
|
||
<div class="outline-text-3" id="text-1-6">
|
||
<p>
|
||
The Cumulative Power Spectrum is the cumulative integral of the Power Spectral Density:
|
||
</p>
|
||
\begin{equation}
|
||
\label{org644886e}
|
||
CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2]
|
||
\end{equation}
|
||
|
||
<p>
|
||
It is also possible to integrate from high frequency to low frequency:
|
||
</p>
|
||
\begin{equation}
|
||
\label{org47d268d}
|
||
CPS_v(f) = \int_f^\infty PSD_v(\nu) d\nu \quad [(m/s)^2]
|
||
\end{equation}
|
||
|
||
<p>
|
||
The Cumulative Power Spectrum taken at frequency \(f\) thus represent the power in the signal in the frequency band \(0\) to \(f\) or \(f\) to \(\infty\) depending on the above definition taken.
|
||
</p>
|
||
|
||
|
||
<p>
|
||
The choice of the integral direction depends on the shape of the PSD.
|
||
If the power is mostly present at low frequencies, it is preferable to use equation \eqref{org47d268d}.
|
||
</p>
|
||
|
||
|
||
<p>
|
||
The Cumulative Amplitude Spectrum is defined as the square root of the Cumulative Power Spectrum:
|
||
\[ CAS_v(f) = \sqrt{CPS_v(f)} = \sqrt{\int_f^\infty PSD_v(\nu) d\nu} \quad [m/s] \]
|
||
</p>
|
||
|
||
<p>
|
||
The Root Mean Square value of the velocity corresponds to the Cumulative Amplitude Spectrum when integrated at all frequencies:
|
||
\[ v_{\text{rms}} = \sqrt{\int_0^\infty PSD_v(\nu) d\nu} = CAS_v(0) \quad [m/s \ \text{rms}] \]
|
||
</p>
|
||
|
||
<p>
|
||
With Matlab, the Cumulative Power Spectrum can be computed with the below formulas and the results are shown in figure <a href="#orgc3a558f">10</a>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> CPS_v = cumtrapz(f, Sv); <span class="org-comment">% Cumulative Power Spectrum from low to high frequencies</span>
|
||
CPS_vv = flip(<span class="org-type">-</span>cumtrapz(flip(f), flip(Sv))); <span class="org-comment">% Cumulative Power Spectrum from high to low frequencies</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> <span class="org-type">figure</span>;
|
||
ax1 = subplot(1, 2, 1);
|
||
hold on;
|
||
plot(f, CPS_v);
|
||
hold off;
|
||
<span class="org-type">set</span>(<span class="org-variable-name">gca</span>, <span class="org-string">'xscale'</span>, <span class="org-string">'log'</span>); <span class="org-type">set</span>(<span class="org-variable-name">gca</span>, <span class="org-string">'yscale'</span>, <span class="org-string">'log'</span>);
|
||
xlabel(<span class="org-string">'Frequency [Hz]'</span>); ylabel(<span class="org-string">'Cumulative Power Spectrum [$(m/s)^2$]'</span>)
|
||
|
||
ax2 = subplot(1, 2, 2);
|
||
hold on;
|
||
plot(f, CPS_vv);
|
||
hold off;
|
||
<span class="org-type">set</span>(<span class="org-variable-name">gca</span>, <span class="org-string">'xscale'</span>, <span class="org-string">'log'</span>); <span class="org-type">set</span>(<span class="org-variable-name">gca</span>, <span class="org-string">'yscale'</span>, <span class="org-string">'log'</span>);
|
||
xlabel(<span class="org-string">'Frequency [Hz]'</span>);
|
||
|
||
linkaxes([ax1,ax2],<span class="org-string">'xy'</span>);
|
||
xlim([0.1, 500]); ylim([1e<span class="org-type">-</span>15, 1e<span class="org-type">-</span>11]);
|
||
</pre>
|
||
</div>
|
||
|
||
|
||
<div id="orgc3a558f" class="figure">
|
||
<p><img src="figs/cps_integral_comp.png" alt="cps_integral_comp.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 10: </span>Cumulative Power Spectrum (<a href="./figs/cps_integral_comp.png">png</a>, <a href="./figs/cps_integral_comp.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org7953d16" class="outline-2">
|
||
<h2 id="org7953d16"><span class="section-number-2">2</span> Time domain signal that approximate a PSD - TF technique</h2>
|
||
<div class="outline-text-2" id="text-2">
|
||
<p>
|
||
<a id="org03ec73d"></a>
|
||
</p>
|
||
</div>
|
||
|
||
<div id="outline-container-org4a1b3f8" class="outline-3">
|
||
<h3 id="org4a1b3f8"><span class="section-number-3">2.1</span> Signal’s PSD</h3>
|
||
<div class="outline-text-3" id="text-2-1">
|
||
<p>
|
||
We load the PSD of the signal we wish to replicate.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> load(<span class="org-string">'./mat/dist_psd.mat'</span>, <span class="org-string">'dist_f'</span>);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We remove the first value with very high PSD.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> dist_f.f = dist_f.f(3<span class="org-type">:</span>end);
|
||
dist_f.psd_gm = dist_f.psd_gm(3<span class="org-type">:</span>end);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The PSD of the signal is shown on figure <a href="#fig:psd_ground_motion">fig:psd_ground_motion</a>.
|
||
</p>
|
||
|
||
<div id="org1b6f6a4" class="figure">
|
||
<p><img src="figs/psd_ground_motion.png" alt="psd_ground_motion.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 11: </span>PSD of the signal (<a href="./figs/psd_ground_motion.png">png</a>, <a href="./figs/psd_ground_motion.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgeba7bdf" class="outline-3">
|
||
<h3 id="orgeba7bdf"><span class="section-number-3">2.2</span> Transfer Function that approximate the ASD</h3>
|
||
<div class="outline-text-3" id="text-2-2">
|
||
<p>
|
||
Using <code>sisotool</code> or any other tool, we create a transfer function \(G\) such that its magnitude is close to the Amplitude Spectral Density \(\Gamma_x = \sqrt{S_x}\):
|
||
\[ |G(j\omega)| \approx \Gamma_x(\omega) \]
|
||
</p>
|
||
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> G_gm = 0.002<span class="org-type">*</span>(s<span class="org-type">^</span>2 <span class="org-type">+</span> 3.169<span class="org-type">*</span>s <span class="org-type">+</span> 27.74)<span class="org-type">/</span>(s<span class="org-type">*</span>(s<span class="org-type">+</span>32.73)<span class="org-type">*</span>(s<span class="org-type">+</span>8.829)<span class="org-type">*</span>(s<span class="org-type">+</span>7.983)<span class="org-type">^</span>2);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We compare the ASD \(\Gamma_x(\omega)\) and the magnitude of the generated transfer function \(|G(j\omega)|\) in figure <a href="#org4f27987">12</a>.
|
||
</p>
|
||
|
||
<div id="org4f27987" class="figure">
|
||
<p><img src="figs/asd_and_tf_compare.png" alt="asd_and_tf_compare.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 12: </span>Comparison of the ASD and of the transfer function’s magnitude (<a href="./figs/asd_and_tf_compare.png">png</a>, <a href="./figs/asd_and_tf_compare.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgbcffeed" class="outline-3">
|
||
<h3 id="orgbcffeed"><span class="section-number-3">2.3</span> Generated Time domain signal</h3>
|
||
<div class="outline-text-3" id="text-2-3">
|
||
<p>
|
||
We know that a signal \(u\) going through a LTI system \(G\) (figure <a href="#org88c197a">13</a>) will have its ASD modified according to the following equation:
|
||
</p>
|
||
\begin{equation}
|
||
\Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{uu}(\omega)
|
||
\end{equation}
|
||
|
||
|
||
<div id="org88c197a" class="figure">
|
||
<p><img src="figs/velocity_to_voltage_psd.png" alt="velocity_to_voltage_psd.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 13: </span>Schematic of the instrumentation used for the measurement</p>
|
||
</div>
|
||
|
||
<p>
|
||
Thus, if we create a random signal with an ASD equal to one at all frequency and we pass this signal through the previously defined LTI transfer function <code>G_gm</code>, we should obtain a signal with an ASD that approximate the original ASD.
|
||
</p>
|
||
|
||
<p>
|
||
To obtain a random signal with an ASD equal to one, we use the following code.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> Fs = 2<span class="org-type">*</span>dist_f.f(end); <span class="org-comment">% Sampling Frequency [Hz]</span>
|
||
Ts = 1<span class="org-type">/</span>Fs; <span class="org-comment">% Sampling Time [s]</span>
|
||
|
||
t = 0<span class="org-type">:</span>Ts<span class="org-type">:</span>500; <span class="org-comment">% Time Vector [s]</span>
|
||
u = sqrt(Fs<span class="org-type">/</span>2)<span class="org-type">*</span>randn(length(t), 1); <span class="org-comment">% Signal with an ASD equal to one</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We then use <code>lsim</code> to compute \(y\) as shown in figure <a href="#org88c197a">13</a>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> y = lsim(G_gm, u, t);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The obtained time domain signal is shown in figure <a href="#org7f4ff26">14</a>.
|
||
</p>
|
||
|
||
<div id="org7f4ff26" class="figure">
|
||
<p><img src="figs/time_domain_u.png" alt="time_domain_u.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 14: </span>Obtained time domain signal \(y(t)\) (<a href="./figs/time_domain_u.png">png</a>, <a href="./figs/time_domain_u.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org676e750" class="outline-3">
|
||
<h3 id="org676e750"><span class="section-number-3">2.4</span> Comparison of the Power Spectral Densities</h3>
|
||
<div class="outline-text-3" id="text-2-4">
|
||
<p>
|
||
We now compute the Power Spectral Density of the computed time domain signal \(y\).
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(y);
|
||
na = 16;
|
||
win = hanning(floor(nx<span class="org-type">/</span>na));
|
||
|
||
[pxx, f] = pwelch(y, win, 0, [], Fs);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Finally, we compare the PSD of the original signal and the obtained signal on figure <a href="#fig:psd_comparison">fig:psd_comparison</a>.
|
||
</p>
|
||
|
||
<div id="orgaa3a04e" class="figure">
|
||
<p><img src="figs/compare_psd_tf_technique.png" alt="compare_psd_tf_technique.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 15: </span>Comparison of the original PSD and the PSD of the computed time domain signal (<a href="./figs/compare_psd_tf_technique.png">png</a>, <a href="./figs/compare_psd_tf_technique.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org8669bb0" class="outline-3">
|
||
<h3 id="org8669bb0"><span class="section-number-3">2.5</span> Simulink</h3>
|
||
<div class="outline-text-3" id="text-2-5">
|
||
<p>
|
||
One advantage of this technique is that it can be easily integrated into simulink.
|
||
</p>
|
||
|
||
<p>
|
||
The corresponding schematic is shown in figure <a href="#orgcee83cb">16</a> where the block <code>Band-Limited White Noise</code> is used to generate a random signal with a PSD equal to one (parameter <code>Noise Power</code> is set to 1).
|
||
</p>
|
||
|
||
<p>
|
||
Then, the signal generated pass through the transfer function representing the wanted ASD.
|
||
</p>
|
||
|
||
|
||
<div id="orgcee83cb" class="figure">
|
||
<p><img src="figs/simulink_psd_generate.png" alt="simulink_psd_generate.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 16: </span>Simulink Schematic</p>
|
||
</div>
|
||
|
||
<p>
|
||
We simulate the system shown in figure <a href="#orgcee83cb">16</a>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> out = <span class="org-matlab-simulink-keyword">sim</span>(<span class="org-string">'matlab/generate_signal_psd.slx'</span>);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
And we compute the PSD of the generated signal.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(out.u_gm.Data);
|
||
na = 8;
|
||
win = hanning(floor(nx<span class="org-type">/</span>na));
|
||
|
||
[pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Finally, we compare the PSD of the generated signal with the original PSD in figure <a href="#org6297d46">17</a>.
|
||
</p>
|
||
|
||
<div id="org6297d46" class="figure">
|
||
<p><img src="figs/compare_psd_original_simulink.png" alt="compare_psd_original_simulink.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 17: </span>Comparison of the obtained signal’s PSD and original PSD (<a href="./figs/compare_psd_original_simulink.png">png</a>, <a href="./figs/compare_psd_original_simulink.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org7073f62" class="outline-2">
|
||
<h2 id="org7073f62"><span class="section-number-2">3</span> Time domain signal that approximate a PSD - IFFT technique</h2>
|
||
<div class="outline-text-2" id="text-3">
|
||
<p>
|
||
<a id="org401f4ee"></a>
|
||
</p>
|
||
<p>
|
||
The technique comes from (<a href="#citeproc_bib_item_1">Preumont 1994</a>) (section 12.11).
|
||
It is used to compute a periodic signal that has any Power Spectral Density defined.
|
||
It makes used of the Unversed Fast Fourier Transform (IFFT).
|
||
</p>
|
||
</div>
|
||
|
||
<div id="outline-container-orge49d520" class="outline-3">
|
||
<h3 id="orge49d520"><span class="section-number-3">3.1</span> Signal’s PSD</h3>
|
||
<div class="outline-text-3" id="text-3-1">
|
||
<p>
|
||
We load the PSD of the signal we wish to replicate.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> load(<span class="org-string">'./mat/dist_psd.mat'</span>, <span class="org-string">'dist_f'</span>);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We remove the first value with very high PSD.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> dist_f.f = dist_f.f(3<span class="org-type">:</span>end);
|
||
dist_f.psd_gm = dist_f.psd_gm(3<span class="org-type">:</span>end);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The PSD of the signal is shown on figure <a href="#fig:psd_original">fig:psd_original</a>.
|
||
</p>
|
||
|
||
<div id="org3dfcdb9" class="figure">
|
||
<p><img src="figs/psd_original.png" alt="psd_original.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 18: </span>PSD of the original signal (<a href="./figs/psd_original.png">png</a>, <a href="./figs/psd_original.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgcb7dbc6" class="outline-3">
|
||
<h3 id="orgcb7dbc6"><span class="section-number-3">3.2</span> Algorithm</h3>
|
||
<div class="outline-text-3" id="text-3-2">
|
||
<p>
|
||
We define some parameters that will be used in the algorithm.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> Fs = 2<span class="org-type">*</span>dist_f.f(end); <span class="org-comment">% Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]</span>
|
||
N = 2<span class="org-type">*</span>length(dist_f.f); <span class="org-comment">% Number of Samples match the one of the wanted PSD</span>
|
||
T0 = N<span class="org-type">/</span>Fs; <span class="org-comment">% Signal Duration [s]</span>
|
||
df = 1<span class="org-type">/</span>T0; <span class="org-comment">% Frequency resolution of the DFT [Hz]</span>
|
||
<span class="org-comment">% Also equal to (dist_f.f(2)-dist_f.f(1))</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We then specify the wanted PSD.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> phi = dist_f.psd_gm;
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We create amplitudes corresponding to wanted PSD.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> C = zeros(N<span class="org-type">/</span>2,1);
|
||
<span class="org-keyword">for</span> <span class="org-variable-name"><span class="org-constant">i</span></span> = <span class="org-constant">1:N/2</span>
|
||
C(<span class="org-constant">i</span>) = sqrt(phi(<span class="org-constant">i</span>)<span class="org-type">*</span>df);
|
||
<span class="org-keyword">end</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Finally, we add some random phase to <code>C</code>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> theta = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>rand(N<span class="org-type">/</span>2,1); <span class="org-comment">% Generate random phase [rad]</span>
|
||
|
||
Cx = [0 ; C<span class="org-type">.*</span>complex(cos(theta),sin(theta))];
|
||
Cx = [Cx; flipud(conj(Cx(2<span class="org-type">:</span>end)))];;
|
||
</pre>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgf14004b" class="outline-3">
|
||
<h3 id="orgf14004b"><span class="section-number-3">3.3</span> Obtained Time Domain Signal</h3>
|
||
<div class="outline-text-3" id="text-3-3">
|
||
<p>
|
||
The time domain data is generated by an inverse FFT.
|
||
</p>
|
||
|
||
<p>
|
||
The <code>ifft</code> Matlab does not take into account the sampling frequency, thus we need to normalize the signal.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> u = N<span class="org-type">/</span>sqrt(2)<span class="org-type">*</span>ifft(Cx); <span class="org-comment">% Normalisation of the IFFT</span>
|
||
t = linspace(0, T0, N<span class="org-type">+</span>1); <span class="org-comment">% Time Vector [s]</span>
|
||
</pre>
|
||
</div>
|
||
|
||
|
||
<div id="orga1f1d0c" class="figure">
|
||
<p><img src="figs/signal_time_domain.png" alt="signal_time_domain.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 19: </span>Obtained signal in the time domain (<a href="./figs/signal_time_domain.png">png</a>, <a href="./figs/signal_time_domain.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org5d770d3" class="outline-3">
|
||
<h3 id="org5d770d3"><span class="section-number-3">3.4</span> PSD Comparison</h3>
|
||
<div class="outline-text-3" id="text-3-4">
|
||
<p>
|
||
We duplicate the time domain signal to have a longer signal and thus a more precise PSD result.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> u_rep = repmat(u, 10, 1);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We compute the PSD of the obtained signal with the following commands.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(u_rep);
|
||
na = 16;
|
||
win = hanning(floor(nx<span class="org-type">/</span>na));
|
||
|
||
[pxx, f] = pwelch(u_rep, win, 0, [], Fs);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Finally, we compare the PSD of the original signal and the obtained signal on figure <a href="#fig:psd_comparison">fig:psd_comparison</a>.
|
||
</p>
|
||
|
||
<div id="org8eed0cb" class="figure">
|
||
<p><img src="figs/psd_comparison.png" alt="psd_comparison.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 20: </span>Comparison of the PSD of the original signal and the PSD of the obtained signal (<a href="./figs/psd_comparison.png">png</a>, <a href="./figs/psd_comparison.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org9515c46" class="outline-2">
|
||
<h2 id="org9515c46"><span class="section-number-2">4</span> Compute the Noise level and Signal level from PSD</h2>
|
||
<div class="outline-text-2" id="text-4">
|
||
<p>
|
||
<a id="orga227890"></a>
|
||
</p>
|
||
<p>
|
||
We here make use of the Power Spectral Density to estimate either the noise level or the amplitude of a deterministic signal.
|
||
Everything is explained in (<a href="#citeproc_bib_item_2">Schmid 2012</a>) sections 5 and 6.
|
||
</p>
|
||
</div>
|
||
|
||
<div id="outline-container-org073b90d" class="outline-3">
|
||
<h3 id="org073b90d"><span class="section-number-3">4.1</span> Time Domain Signal</h3>
|
||
<div class="outline-text-3" id="text-4-1">
|
||
<p>
|
||
Let’s first define the number of sample and the sampling time.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> N = 10000; <span class="org-comment">% Number of Sample</span>
|
||
dt = 0.001; <span class="org-comment">% Sampling Time [s]</span>
|
||
|
||
t = dt<span class="org-type">*</span>(0<span class="org-type">:</span>1<span class="org-type">:</span>N<span class="org-type">-</span>1)<span class="org-type">'</span>; <span class="org-comment">% Time vector [s]</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We generate of signal that consist of:
|
||
</p>
|
||
<ul class="org-ul">
|
||
<li>a white noise with an RMS value equal to <code>anoi</code></li>
|
||
<li>two sinusoidal signals</li>
|
||
</ul>
|
||
|
||
<p>
|
||
The parameters are defined below.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> asig = 0.1; <span class="org-comment">% Amplitude of the signal [V]</span>
|
||
fsig = 10; <span class="org-comment">% Frequency of the signal [Hz]</span>
|
||
|
||
ahar = 0.5; <span class="org-comment">% Amplitude of the harmonic [V]</span>
|
||
fhar = 50; <span class="org-comment">% Frequency of the harmonic [Hz]</span>
|
||
|
||
anoi = 1e<span class="org-type">-</span>3; <span class="org-comment">% RMS value of the noise</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The signal \(x\) is generated with the following code and is shown in figure <a href="#org284b762">21</a>.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> x = anoi<span class="org-type">*</span>randn(N, 1) <span class="org-type">+</span> asig<span class="org-type">*</span>sin((2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>fsig)<span class="org-type">*</span>t) <span class="org-type">+</span> ahar<span class="org-type">*</span>sin((2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>fhar)<span class="org-type">*</span>t);
|
||
</pre>
|
||
</div>
|
||
|
||
|
||
<div id="org284b762" class="figure">
|
||
<p><img src="figs/time_domain_x_zoom.png" alt="time_domain_x_zoom.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 21: </span>Time Domain Signal (<a href="./figs/time_domain_x_zoom.png">png</a>, <a href="./figs/time_domain_x_zoom.pdf">pdf</a>)</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org86ba021" class="outline-3">
|
||
<h3 id="org86ba021"><span class="section-number-3">4.2</span> Estimation of the magnitude of a deterministic signal</h3>
|
||
<div class="outline-text-3" id="text-4-2">
|
||
<p>
|
||
Let’s compute the PSD of the signal using the <code>blackmanharris</code> window.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nx = length(x);
|
||
na = 8;
|
||
win = blackmanharris(floor(nx<span class="org-type">/</span>na));
|
||
|
||
[pxx, f] = pwelch(x, win, 0, [], 1<span class="org-type">/</span>dt);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
Normalization of the PSD.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> CG = sum(win)<span class="org-type">/</span>(nx<span class="org-type">/</span>na);
|
||
NG = sum(win<span class="org-type">.^</span>2)<span class="org-type">/</span>(nx<span class="org-type">/</span>na);
|
||
fbin = f(2) <span class="org-type">-</span> f(1);
|
||
|
||
pxx_norm = pxx<span class="org-type">*</span>(NG<span class="org-type">*</span>fbin<span class="org-type">/</span>CG<span class="org-type">^</span>2);
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We determine the frequency bins corresponding to the frequency of the signals.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> isig = round(fsig<span class="org-type">/</span>fbin)<span class="org-type">+</span>1;
|
||
ihar = round(fhar<span class="org-type">/</span>fbin)<span class="org-type">+</span>1;
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
The theoretical RMS value of the signal is:
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> srmt = asig<span class="org-type">/</span>sqrt(2); <span class="org-comment">% Theoretical value of signal magnitude</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
And we estimate the RMS value of the signal by either integrating the PSD around the frequency of the signal or by just taking the maximum value.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> srms = sqrt(sum(pxx(isig<span class="org-type">-</span>5<span class="org-type">:</span>isig<span class="org-type">+</span>5)<span class="org-type">*</span>fbin)); <span class="org-comment">% Signal spectrum integrated</span>
|
||
srmsp = sqrt(pxx_norm(isig) <span class="org-type">*</span> NG<span class="org-type">*</span>fbin<span class="org-type">/</span>CG<span class="org-type">^</span>2); <span class="org-comment">% Maximum read off spectrum</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<table border="2" cellspacing="0" cellpadding="6" rules="groups" frame="hsides">
|
||
|
||
|
||
<colgroup>
|
||
<col class="org-left" />
|
||
|
||
<col class="org-right" />
|
||
</colgroup>
|
||
<thead>
|
||
<tr>
|
||
<th scope="col" class="org-left"> </th>
|
||
<th scope="col" class="org-right">Signal Magnitude [rms]</th>
|
||
</tr>
|
||
</thead>
|
||
<tbody>
|
||
<tr>
|
||
<td class="org-left">Theoretical</td>
|
||
<td class="org-right">0.0707</td>
|
||
</tr>
|
||
|
||
<tr>
|
||
<td class="org-left">Integrated</td>
|
||
<td class="org-right">0.0707</td>
|
||
</tr>
|
||
|
||
<tr>
|
||
<td class="org-left">Maximum</td>
|
||
<td class="org-right">0.0529</td>
|
||
</tr>
|
||
</tbody>
|
||
</table>
|
||
|
||
<p>
|
||
We see that the integrated value gives a good approximate of the true RMS value whereas the maximum value give an poor approximate.
|
||
This effect is called <i>scallop loss</i>.
|
||
Thus, always the integrated method should be used.
|
||
</p>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-org833e18d" class="outline-3">
|
||
<h3 id="org833e18d"><span class="section-number-3">4.3</span> Estimation of the noise level</h3>
|
||
<div class="outline-text-3" id="text-4-3">
|
||
<p>
|
||
The noise level can also be computed using the integration method.
|
||
</p>
|
||
|
||
<p>
|
||
The theoretical RMS noise value is.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> nth = anoi<span class="org-type">/</span>sqrt(max(f)) <span class="org-comment">% Theoretical value [V/sqrt(Hz)]</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<p>
|
||
We can estimate this RMS value by integrating the PSD at frequencies where the power of the noise signal is above the power of the other signals.
|
||
</p>
|
||
<div class="org-src-container">
|
||
<pre class="src src-matlab"> navg = sqrt(mean(pxx_norm([ihar<span class="org-type">+</span>10<span class="org-type">:</span>end]))) <span class="org-comment">% pwelch output averaged</span>
|
||
</pre>
|
||
</div>
|
||
|
||
<table border="2" cellspacing="0" cellpadding="6" rules="groups" frame="hsides">
|
||
|
||
|
||
<colgroup>
|
||
<col class="org-left" />
|
||
|
||
<col class="org-right" />
|
||
</colgroup>
|
||
<thead>
|
||
<tr>
|
||
<th scope="col" class="org-left"> </th>
|
||
<th scope="col" class="org-right">test</th>
|
||
</tr>
|
||
</thead>
|
||
<tbody>
|
||
<tr>
|
||
<td class="org-left">Theoretical</td>
|
||
<td class="org-right">4.472e-05</td>
|
||
</tr>
|
||
|
||
<tr>
|
||
<td class="org-left">Average</td>
|
||
<td class="org-right">4.337e-05</td>
|
||
</tr>
|
||
</tbody>
|
||
</table>
|
||
|
||
<p>
|
||
The estimate of the noise level is quite good.
|
||
</p>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div id="outline-container-orgd325fcb" class="outline-2">
|
||
<h2 id="orgd325fcb"><span class="section-number-2">5</span> Further Notes</h2>
|
||
<div class="outline-text-2" id="text-5">
|
||
</div>
|
||
<div id="outline-container-orgf74d5dc" class="outline-3">
|
||
<h3 id="orgf74d5dc"><span class="section-number-3">5.1</span> PSD of ADC quantization noise</h3>
|
||
<div class="outline-text-3" id="text-5-1">
|
||
<p>
|
||
This is taken from <a href="https://www.allaboutcircuits.com/technical-articles/quantization-nois-amplitude-quantization-error-analog-to-digital-converters/">here</a>.
|
||
</p>
|
||
|
||
<p>
|
||
Let’s note:
|
||
</p>
|
||
<ul class="org-ul">
|
||
<li>\(q\) is the corresponding value in [V] of the least significant bit (LSB)</li>
|
||
<li>\(\Delta V\) is the full range of the ADC in [V]</li>
|
||
<li>\(n\) is the number of ADC’s bits</li>
|
||
<li>\(f_s\) is the sample frequency in [Hz]</li>
|
||
</ul>
|
||
|
||
<p>
|
||
Let’s suppose that the ADC is ideal.
|
||
The only noise comes from the quantization error.
|
||
Interestingly, the noise amplitude is uniformly distributed.
|
||
</p>
|
||
|
||
<p>
|
||
The quantization noise can take a value between \(\pm q/2\), and the probability density function is constant in this range (i.e., it’s a uniform distribution).
|
||
Since the integral of the probability density function is equal to one, its value will be \(1/q\) for \(-q/2 < e < q/2\) (Fig. <a href="#orge2c7d6f">22</a>).
|
||
</p>
|
||
|
||
|
||
<div id="orge2c7d6f" class="figure">
|
||
<p><img src="figs/probability_density_function_adc.png" alt="probability_density_function_adc.png" />
|
||
</p>
|
||
<p><span class="figure-number">Figure 22: </span>Probability density function \(p(e)\) of the ADC error \(e\)</p>
|
||
</div>
|
||
|
||
<p>
|
||
Now, we can calculate the time average power of the quantization noise as
|
||
</p>
|
||
\begin{equation}
|
||
P_q = \int_{-q/2}^{q/2} e^2 p(e) de = \frac{q^2}{12}
|
||
\end{equation}
|
||
|
||
|
||
<p>
|
||
The other important parameter of a noise source is the power spectral density (PSD), which indicates how the noise power spreads in different frequency bands.
|
||
To find the power spectral density, we need to calculate the Fourier transform of the autocorrelation function of the noise.
|
||
</p>
|
||
|
||
<p>
|
||
Assuming that the noise samples are not correlated with one another, we can approximate the autocorrelation function with a delta function in the time domain.
|
||
Since the Fourier transform of a delta function is equal to one, the <b>power spectral density will be frequency independent</b>.
|
||
Therefore, the quantization noise is white noise with total power equal to \(P_q = \frac{q^2}{12}\).
|
||
</p>
|
||
|
||
<p>
|
||
Thus, the two-sided PSD (from \(\frac{-f_s}{2}\) to \(\frac{f_s}{2}\)), we should divide the noise power \(P_q\) by \(f_s\):
|
||
</p>
|
||
|
||
\begin{equation}
|
||
\int_{-f_s/2}^{f_s/2} \Gamma(f) d f = f_s \Gamma = \frac{q^2}{12}
|
||
\end{equation}
|
||
|
||
<p>
|
||
Finally:
|
||
</p>
|
||
\begin{equation}
|
||
\begin{aligned}
|
||
\Gamma &= \frac{q^2}{12 f_s} \\
|
||
&= \frac{\left(\frac{\Delta V}{2^n}\right)^2}{12 f_s} \text{ in } \left[ \frac{V^2}{Hz} \right]
|
||
\end{aligned}
|
||
\end{equation}
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<p>
|
||
|
||
</p>
|
||
|
||
<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><h2 class='citeproc-org-bib-h2'>Bibliography</h2>
|
||
<div class="csl-bib-body">
|
||
<div class="csl-entry"><a name="citeproc_bib_item_1"></a>Preumont, André. 1994. <i>Random Vibration and Spectral Analysis</i>. Solid Mechanics and Its Applications. Springer Netherlands. <a href="https://doi.org/10.1007/978-94-017-2840-9">https://doi.org/10.1007/978-94-017-2840-9</a>.</div>
|
||
<div class="csl-entry"><a name="citeproc_bib_item_2"></a>Schmid, Hanspeter. 2012. “How to Use the FFT and Matlab’s Pwelch Function for Signal and Noise Simulations and Measurements.” <i>Institute of Microelectronics</i>.</div>
|
||
</div>
|
||
</div>
|
||
<div id="postamble" class="status">
|
||
<p class="author">Author: Dehaeze Thomas</p>
|
||
<p class="date">Created: 2021-04-26 lun. 17:04</p>
|
||
</div>
|
||
</body>
|
||
</html>
|