Compute Spectral Densities of signals with Matlab
Table of Contents
-
-
- 1. Sensitivity of the instrumentation -
- 2. Convert the time domain from volts to velocity -
- 3. Power Spectral Density and Amplitude Spectral Density -
- 4. Modification of a signal's Power Spectral Density when going through an LTI system -
- 5. From PSD of the velocity to the PSD of the displacement +
- 1. Spectral Analysis - Basics
+
-
+
- 1.1. PSD of ADC quantization noise +
- 1.2. Sensitivity of the instrumentation +
- 1.3. Convert the time domain from volts to velocity +
- 1.4. Power Spectral Density and Amplitude Spectral Density +
- 1.5. Modification of a signal's Power Spectral Density when going through an LTI system +
- 1.6. From PSD of the velocity to the PSD of the displacement +
- 1.7. Cumulative Power/Amplitude Spectrum +
- 1.8. TODO Add best practices from the Article and simple snippet that works +
+ - 2. Technique 1 : Approximation with a transfer function + + +
- 3. Technique 2 : IFFT + + +
- 4. TODO Compute the Noise level and Signal level from PSD + +
-This document presents the mathematics as well as the matlab scripts to do the spectral analysis of a measured signal. +This document presents the mathematics as well as the matlab scripts to do various spectral analysis on a measured signal.
++Some matlab documentation about Spectral Analysis can be found here. +
+ ++First, in section 1, some basics of spectral analysis are presented. +
+ ++In some cases, we want to generate a time domain signal with defined Power Spectral Density. +Two methods are presented in sections 2 and 3. +
+ ++Finally, some notes are done on how to compute the noise level and signal level from a given Power Spectral Density in section 4. +
+ +
+
1 Spectral Analysis - Basics
+
+
-Typically this signal is coming from an inertial sensor, a force sensor or any other sensor.
@@ -300,10 +348,85 @@ Typically this signal is coming from an inertial sensor, a force sensor or any oWe here take the example of a signal coming from a Geophone measurement the vertical velocity of the floor at the ESRF.
+
-
1 Sensitivity of the instrumentation
-
+
+
+
+1.1 PSD of ADC quantization noise
+
+
++This is taken from here. +
+ ++Let's note: +
+-
+
- \(q\) is the corresponding value in [V] of the least significant bit (LSB) +
- \(\Delta V\) is the full range of the ADC in [V] +
- \(n\) is the number of ADC's bits +
- \(f_s\) is the sample frequency in [Hz] +
+Let's suppose that the ADC is ideal. +The only noise comes from the quantization error. +Interestingly, the noise amplitude is uniformly distributed. +
+ ++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. 1). +
+ + +
+
+
+
+
Figure 1: Probability density function \(p(e)\) of the ADC error \(e\)
++Now, we can calculate the time average power of the quantization noise as +
+\begin{equation} + P_q = \int_{-q/2}^{q/2} e^2 p(e) de = \frac{q^2}{12} +\end{equation} + + ++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. +
+ ++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 power spectral density will be frequency independent. +Therefore, the quantization noise is white noise with total power equal to \(P_q = \frac{q^2}{12}\). +
+ ++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\): +
+ +\begin{equation} + \int_{-f_s/2}^{f_s/2} \Gamma(f) d f = f_s \Gamma = \frac{q^2}{12} +\end{equation} + ++Finally: +
+\begin{equation} + \begin{align} + \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{align} +\end{equation} +
+
-
1.2 Sensitivity of the instrumentation
+The measured signal \(x\) by the ADC is in Volts. The corresponding real velocity \(v\) in m/s. @@ -314,17 +437,17 @@ To obtain the real quantity as measured by the sensor, one have to know the sens
-
+
Figure 1: Schematic of the instrumentation used for the measurement
+Figure 2: Schematic of the instrumentation used for the measurement
-
2 Convert the time domain from volts to velocity
-
+
+
1.3 Convert the time domain from volts to velocity
+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] \] @@ -348,11 +471,11 @@ And the gain of the amplifier is 1000: \(G_m(s) = 1000\).
-If \({G_m(s)}^{-1} {G_g(s)}^{-1}\) is proper, we can simulate this dynamical system to go from the voltage to the velocity units (figure 2). +If \({G_m(s)}^{-1} {G_g(s)}^{-1}\) is proper, we can simulate this dynamical system to go from the voltage to the velocity units (figure 3).
-
-data = load('mat/data_028.mat', 'data'); data = data.data; +data = load('mat/data_028.mat', 'data'); data = data.data; t = data(:, 3); % [s] x = data(:, 1)-mean(data(:, 1)); % The offset if removed (coming from the voltage amplifier) [v] @@ -362,17 +485,17 @@ dt = t( +-
Figure 2: Schematic of the instrumentation used for the measurement
+Figure 3: Schematic of the instrumentation used for the measurement
We simulate this system with matlab:
-@@ -382,22 +505,22 @@ And we plot the obtained velocityv = lsim(inv(Gg*Gm), v, t); +v = lsim(inv(Gg*Gm), x, t);-figure; plot(t, v); -xlabel("Time [s]"); ylabel("Velocity [m/s]"); +xlabel("Time [s]"); ylabel("Velocity [m/s]");+-
Figure 3: Measured Velocity
+Figure 4: Measured Velocity
-
3 Power Spectral Density and Amplitude Spectral Density
-
+
+
-1.4 Power Spectral Density and Amplitude Spectral Density
+We now have the velocity in the time domain: \[ v(t)\ [m/s] \] @@ -437,27 +560,27 @@ ylabel('Power Spectral Densi The Amplitude Spectral Density (ASD) is the square root of the Power Spectral Density:
\begin{equation} - \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right] + \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right] \end{equation}figure; loglog(f, sqrt(Sv)); xlabel('Frequency [Hz]'); -ylabel('Amplitude Spectral Density $\left[\frac{m/s}{\sqrt{Hz}}\right]$') +ylabel('Amplitude Spectral Density $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
-
4 Modification of a signal's Power Spectral Density when going through an LTI system
-
+
+
1.5 Modification of a signal's Power Spectral Density when going through an LTI system
+
-
-
+
Figure 4: Schematic of the instrumentation used for the measurement
+Figure 5: Schematic of the instrumentation used for the measurement
@@ -476,14 +599,14 @@ And we also have:
->
+
+** Introduction :ignore:
Typically this signal is coming from an inertial sensor, a force sensor or any other sensor.
We here take the example of a signal coming from a Geophone measurement the vertical velocity of the floor at the ESRF.
-* Matlab Init :noexport:ignore:
+** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<>
#+end_src
@@ -49,7 +63,68 @@ We here take the example of a signal coming from a Geophone measurement the vert
<>
#+end_src
-* Sensitivity of the instrumentation
+** PSD of ADC quantization noise
+This is taken from [[https://www.allaboutcircuits.com/technical-articles/quantization-nois-amplitude-quantization-error-analog-to-digital-converters/][here]].
+
+Let's note:
+- $q$ is the corresponding value in [V] of the least significant bit (LSB)
+- $\Delta V$ is the full range of the ADC in [V]
+- $n$ is the number of ADC's bits
+- $f_s$ is the sample frequency in [Hz]
+
+Let's suppose that the ADC is ideal.
+The only noise comes from the quantization error.
+Interestingly, the noise amplitude is uniformly distributed.
+
+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. [[fig:probability_density_function_adc]]).
+
+#+begin_src latex :file probability_density_function_adc.pdf :post pdf2svg(file=*this*, ext="png") :exports results
+ \begin{tikzpicture}
+ \path[fill=black!20!white] (-1, 0) |- (1, 1) |- (-1, 0);
+
+ \draw[->] (-2, 0) -- (2, 0) node[above left]{$e$};
+ \draw[->] (0, -0.5) -- (0, 2) node[below right]{$p(e)$};
+
+ \node[below] at (1, 0){$\frac{q}{2}$};
+ \node[below] at (-1, 0){$-\frac{q}{2}$};
+ \node[right] at (1, 1){$\frac{1}{q}$};
+ \end{tikzpicture}
+#+end_src
+
+#+name: fig:probability_density_function_adc
+#+caption: Probability density function $p(e)$ of the ADC error $e$
+#+RESULTS:
+[[file:figs/probability_density_function_adc.png]]
+
+Now, we can calculate the time average power of the quantization noise as
+\begin{equation}
+ P_q = \int_{-q/2}^{q/2} e^2 p(e) de = \frac{q^2}{12}
+\end{equation}
+
+
+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.
+
+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 *power spectral density will be frequency independent*.
+Therefore, the quantization noise is white noise with total power equal to $P_q = \frac{q^2}{12}$.
+
+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$:
+
+\begin{equation}
+ \int_{-f_s/2}^{f_s/2} \Gamma(f) d f = f_s \Gamma = \frac{q^2}{12}
+\end{equation}
+
+Finally:
+\begin{equation}
+ \begin{align}
+ \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{align}
+\end{equation}
+
+** Sensitivity of the instrumentation
The measured signal $x$ by the ADC is in Volts.
The corresponding real velocity $v$ in m/s.
@@ -75,7 +150,7 @@ To obtain the real quantity as measured by the sensor, one have to know the sens
#+RESULTS: fig:velocity_to_voltage
[[file:figs/velocity_to_voltage.png]]
-* Convert the time domain from volts to velocity
+** Convert the time domain from volts to velocity
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] \]
@@ -123,7 +198,7 @@ If ${G_m(s)}^{-1} {G_g(s)}^{-1}$ is proper, we can simulate this dynamical syste
We simulate this system with matlab:
#+begin_src matlab
- v = lsim(inv(Gg*Gm), v, t);
+ v = lsim(inv(Gg*Gm), x, t);
#+end_src
And we plot the obtained velocity
@@ -144,7 +219,7 @@ And we plot the obtained velocity
#+RESULTS: fig:velocity_time
[[file:figs/velocity_time.png]]
-* Power Spectral Density and Amplitude Spectral Density
+** Power Spectral Density and Amplitude Spectral Density
We now have the velocity in the time domain:
\[ v(t)\ [m/s] \]
@@ -181,7 +256,7 @@ We first have to defined a window:
The Amplitude Spectral Density (ASD) is the square root of the Power Spectral Density:
\begin{equation}
- \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right]
+ \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right]
\end{equation}
#+begin_src matlab
@@ -201,7 +276,7 @@ The Amplitude Spectral Density (ASD) is the square root of the Power Spectral De
#+CAPTION: Power Spectral Density of the measured velocity
#+RESULTS: fig:asd_velocity
-* Modification of a signal's Power Spectral Density when going through an LTI system
+** Modification of a signal's Power Spectral Density when going through an LTI system
#+begin_src latex :file velocity_to_voltage_psd.pdf :post pdf2svg(file=*this*, ext="png") :exports results
\begin{tikzpicture}
@@ -227,7 +302,7 @@ And we also have:
\Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{xx}(\omega)
\end{equation}
-* From PSD of the velocity to the PSD of the displacement
+** From PSD of the velocity to the PSD of the displacement
#+begin_src latex :file velocity_to_displacement_psd.pdf :post pdf2svg(file=*this*, ext="png") :exports results
\begin{tikzpicture}
@@ -284,6 +359,7 @@ And we have
Note here that we always have
\[ PSD_x \left(f = \frac{1}{2\pi}\right) = PSD_v \left(f = \frac{1}{2\pi}\right) = PSD_a \left(f = \frac{1}{2\pi}\right), \quad \frac{1}{2\pi} \approx 0.16 [Hz] \]
+** Cumulative Power/Amplitude Spectrum
If we want to compute the Cumulative Power Spectrum:
\[ CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2] \]
@@ -296,8 +372,387 @@ The Cumulative Amplitude Spectrum is then the square root of the Cumulative Powe
Then, we can obtain the Root Mean Square value of the velocity:
\[ v_{\text{rms}} = CAS_v(0) \quad [m/s \ \text{rms}] \]
-#+begin_src matlab
+#+begin_src matlab :results none
+ figure;
+ hold on;
+ plot(f, cumtrapz(f, Sv));
+ hold off;
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+ xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]')
+#+end_src
+In order to integrate from high frequency to low frequency:
+#+begin_src matlab :results none
+ figure;
+ hold on;
+ plot(f, flip(-cumtrapz(flip(f), flip(Sv))));
+ hold off;
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+ xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]')
+#+end_src
+
+** TODO Add best practices from the Article and simple snippet that works
+* Technique 1 : Approximation with a transfer function
+<>
+** Introduction :ignore:
+
+** Matlab Init :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+ <>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+ <>
+#+end_src
+
+** Signal's PSD
+We load the PSD of the signal we wish to replicate.
+#+begin_src matlab
+ load('./mat/dist_psd.mat', 'dist_f');
+#+end_src
+
+We remove the first value with very high PSD.
+#+begin_src matlab
+ dist_f.f = dist_f.f(3:end);
+ dist_f.psd_gm = dist_f.psd_gm(3:end);
+#+end_src
+
+The PSD of the signal is shown on figure ref:fig:psd_ground_motion.
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, dist_f.psd_gm)
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/psd_ground_motion.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+ <>
+#+end_src
+
+#+NAME: fig:psd_ground_motion
+#+CAPTION: PSD of the signal ([[./figs/psd_ground_motion.png][png]], [[./figs/psd_ground_motion.pdf][pdf]])
+[[file:figs/psd_ground_motion.png]]
+
+** Transfer Function that approximate the ASD
+#+begin_src matlab
+ G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2);
+#+end_src
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, sqrt(dist_f.psd_gm))
+ plot(dist_f.f, abs(squeeze(freqresp(G_gm, dist_f.f, 'Hz'))))
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+#+end_src
+
+** Generated Time domain signal
+#+begin_src matlab
+ Fs = 2*dist_f.f(end);
+ Ts = 1/Fs;
+
+ t = 0:Ts:500;
+ u = sqrt(Fs/2)*randn(length(t), 1);
+#+end_src
+
+#+begin_src matlab
+ u_o = lsim(G_gm, u, t);
+#+end_src
+
+#+begin_src matlab
+ figure;
+ plot(t, u_o);
+#+end_src
+
+** Comparison of the Power Spectral Densities
+#+begin_src matlab
+ nx = length(u_o);
+ na = 16;
+ win = hanning(floor(nx/na));
+
+ [pxx, f] = pwelch(u_o, win, 0, [], Fs);
+#+end_src
+
+Finally, we compare the PSD of the original signal and the obtained signal on figure ref:fig:psd_comparison.
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
+ plot(f, pxx, 'DisplayName', 'Computed')
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+ legend('location', 'northeast');
+#+end_src
+
+** Simulink
+
+#+name: fig:simulink_psd_generate
+#+caption: Simulink Schematic
+[[file:figs/simulink_psd_generate.png]]
+
+The parameters for the =Band-Limited White Noise= are:
+- Noise Power: 1
+
+#+begin_src matlab
+ nx = length(out.u_gm.Data);
+ na = 8;
+ win = hanning(floor(nx/na));
+
+ [pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);
+#+end_src
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
+ plot(f, pxx, 'DisplayName', 'Computed')
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+ legend('location', 'northeast');
+#+end_src
+
+* Technique 2 : IFFT
+<>
+** Introduction :ignore:
+The technique comes from cite:preumont94_random_vibrat_spect_analy (section 12.11).
+
+** Matlab Init :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+ <>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+ <>
+#+end_src
+
+** Signal's PSD
+We load the PSD of the signal we wish to replicate.
+#+begin_src matlab
+ load('./mat/dist_psd.mat', 'dist_f');
+#+end_src
+
+We remove the first value with very high PSD.
+#+begin_src matlab
+ dist_f.f = dist_f.f(3:end);
+ dist_f.psd_gm = dist_f.psd_gm(3:end);
+#+end_src
+
+The PSD of the signal is shown on figure ref:fig:psd_original.
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, dist_f.psd_gm)
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/psd_original.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+ <>
+#+end_src
+
+#+NAME: fig:psd_original
+#+CAPTION: PSD of the original signal ([[./figs/psd_original.png][png]], [[./figs/psd_original.pdf][pdf]])
+[[file:figs/psd_original.png]]
+
+** Algorithm
+We define some parameters.
+#+begin_src matlab
+ Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+ N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD
+ T0 = N/Fs; % Signal Duration [s]
+ df = 1/T0; % Frequency resolution of the DFT [Hz]
+ % Also equal to (dist_f.f(2)-dist_f.f(1))
+#+end_src
+
+We then specify the wanted PSD.
+#+begin_src matlab
+ phi = dist_f.psd_gm;
+#+end_src
+
+Create amplitudes corresponding to wanted PSD.
+#+begin_src matlab
+ C = zeros(N/2,1);
+ for i = 1:N/2
+ C(i) = sqrt(phi(i)*df);
+ end
+#+end_src
+
+Add random phase to =C=.
+#+begin_src matlab
+ theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+
+ Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+ Cx = [Cx; flipud(conj(Cx(2:end)))];;
+#+end_src
+
+** Obtained Time Domain Signal
+The time domain data is generated by an inverse FFT.
+We normalize the =ifft= Matlab command.
+#+begin_src matlab
+ u = 1/(sqrt(2)*df*1/Fs)*ifft(Cx); % Normalisation of the IFFT
+ t = linspace(0, T0, N+1); % Time Vector [s]
+#+end_src
+
+#+begin_src matlab :exports none
+ figure;
+ plot(t, u)
+ xlabel('Time [s]');
+ ylabel('Amplitude');
+ xlim([t(1), t(end)]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/signal_time_domain.pdf" :var figsize="normal-wide" :post pdf2svg(file=*this*, ext="png")
+ <>
+#+end_src
+
+#+NAME: fig:signal_time_domain
+#+CAPTION: Obtained signal in the time domain ([[./figs/signal_time_domain.png][png]], [[./figs/signal_time_domain.pdf][pdf]])
+[[file:figs/signal_time_domain.png]]
+
+** PSD Comparison
+We duplicate the time domain signal to have a longer signal and thus a more precise PSD result.
+#+begin_src matlab
+ u_rep = repmat(u, 10, 1);
+#+end_src
+
+We compute the PSD of the obtained signal with the following commands.
+#+begin_src matlab
+ nx = length(u_rep);
+ na = 16;
+ win = hanning(floor(nx/na));
+
+ [pxx, f] = pwelch(u_rep, win, 0, [], Fs);
+#+end_src
+
+Finally, we compare the PSD of the original signal and the obtained signal on figure ref:fig:psd_comparison.
+
+#+begin_src matlab :exports none
+ figure;
+ hold on;
+ plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
+ plot(f, pxx, 'DisplayName', 'Computed')
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+ legend('location', 'northeast');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/psd_comparison.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+ <>
+#+end_src
+
+#+NAME: fig:psd_comparison
+#+CAPTION: Comparison of the PSD of the original signal and the PSD of the obtained signal ([[./figs/psd_comparison.png][png]], [[./figs/psd_comparison.pdf][pdf]])
+[[file:figs/psd_comparison.png]]
+
+* TODO Compute the Noise level and Signal level from PSD
+<>
+** Introduction :ignore:
+
+** Matlab Init :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+ <>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+ <>
+#+end_src
+
+** Computation
+- [ ] Add table to compare the methods
+- [ ] Add some explanations
+
+#+begin_src matlab
+ N = 10000;
+ dt = 0.001;
+
+ t = dt*(0:1:N-1)';
+#+end_src
+
+Parameters of the signal
+#+begin_src matlab
+ asig = 0.8; % Amplitude of the signal [V]
+ fsig = 100; % Frequency of the signal [Hz]
+
+ anoi = 1e-3; % RMS value of the noise
+
+ x = anoi*randn(N, 1) + asig*sin((2*pi*fsig)*t);
+#+end_src
+
+#+begin_src matlab
+ figure;
+ plot(t, x);
+#+end_src
+
+Compute the PSD of the signal.
+#+begin_src matlab
+ nx = length(x);
+ na = 8;
+ win = blackmanharris(floor(nx/na));
+
+ [pxx, f] = pwelch(x, win, 0, [], 1/dt);
+#+end_src
+
+Normalization of the PSD.
+#+begin_src matlab
+ CG = sum(win)/(nx/na);
+ NG = sum(win.^2)/(nx/na);
+ fbin = f(2) - f(1);
+
+ pxx_norm = pxx*(NG*fbin/(CG)^2);
+#+end_src
+
+#+begin_src matlab
+ isig = round(fsig/fbin)+1;
+#+end_src
+
+Estimate the Signal magnitude.
+#+begin_src matlab
+ srmt = asig/sqrt(2) % Theoretical value of signal magnitude
+ srms = sqrt(sum(pxx(isig-5:isig+5)*fbin)) % Signal spectrum integrated
+ srmsp = sqrt(pxx(isig) * NG*fbin/CG^2) % Maximum read off spectrum
+#+end_src
+
+Estimate the noise floor.
+#+begin_src matlab
+ nth = anoi/sqrt(max(f)) % Theoretical value [V/sqrt(Hz)]
+
+ inmax = isig-20;
+ nsum = sqrt(sum(pxx(1:inmax)*fbin)) / sqrt(f(inmax)) % Signal spectrum integrated
+
+ navg = sqrt(mean(pxx(1:inmax))) % pwelch output averaged
+#+end_src
+
+#+begin_src matlab
+ figure;
+ hold on;
+ plot(f, pxx)
+ plot(f, pxx_norm)
+ hold off;
+ xlabel('Frequency [Hz]');
+ ylabel('Power Spectral Density');
+ set(gca, 'xscale', 'log');
+ set(gca, 'yscale', 'log');
#+end_src
* Bibliography :ignore:
diff --git a/mat/dist_psd.mat b/mat/dist_psd.mat
new file mode 100644
index 0000000..a82319c
Binary files /dev/null and b/mat/dist_psd.mat differ
diff --git a/matlab/generate_signal_psd.slx b/matlab/generate_signal_psd.slx
new file mode 100644
index 0000000..8c9dbe1
Binary files /dev/null and b/matlab/generate_signal_psd.slx differ
diff --git a/ref.bib b/ref.bib
index 243d8e1..00e9771 100644
--- a/ref.bib
+++ b/ref.bib
@@ -1,8 +1,20 @@
@article{schmid12_how_to_use_fft_matlab,
author = {Schmid, Hanspeter},
- title = {How To Use the Fft and Matlab's Pwelch Function for Signal and
+ title = {How To Use the FFT and Matlab's Pwelch Function for Signal and
Noise Simulations and Measurements},
journal = {Institute of Microelectronics},
year = 2012,
keywords = {},
}
+
+@book{preumont94_random_vibrat_spect_analy,
+ author = {Andr{\'e} Preumont},
+ title = {Random Vibration and Spectral Analysis},
+ year = 1994,
+ publisher = {Springer Netherlands},
+ url = {https://doi.org/10.1007/978-94-017-2840-9},
+ DATE_ADDED = {Mon Dec 2 08:52:51 2019},
+ doi = {10.1007/978-94-017-2840-9},
+ pages = {nil},
+ series = {Solid Mechanics and Its Applications},
+}
\ No newline at end of file
5 From PSD of the velocity to the PSD of the displacement
-
+
diff --git a/index.org b/index.org
index 4264ef2..798b576 100644
--- a/index.org
+++ b/index.org
@@ -12,7 +12,7 @@
#+HTML_HEAD:
#+HTML_HEAD:
-#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
+#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/thesis/latex/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
@@ -34,13 +34,27 @@
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:
-This document presents the mathematics as well as the matlab scripts to do the spectral analysis of a measured signal.
+* Introduction :ignore:
+This document presents the mathematics as well as the matlab scripts to do various spectral analysis on a measured signal.
+Some matlab documentation about Spectral Analysis can be found [[https://fr.mathworks.com/help/signal/ug/spectral-analysis.html][here]].
+
+First, in section [[sec:spectral_analysis_basics]], some basics of spectral analysis are presented.
+
+In some cases, we want to generate a time domain signal with defined Power Spectral Density.
+Two methods are presented in sections [[sec:approximate_tf]] and [[sec:approximate_ifft]].
+
+Finally, some notes are done on how to compute the noise level and signal level from a given Power Spectral Density in section [[sec:compute_psd_levels]].
+
+* Spectral Analysis - Basics
+<
+
+1.6 From PSD of the velocity to the PSD of the displacement
+
-
+
+
+
Figure 5: Schematic of the instrumentation used for the measurement
+Figure 6: Schematic of the instrumentation used for the measurement
@@ -545,7 +668,12 @@ And we have Note here that we always have \[ PSD_x \left(f = \frac{1}{2\pi}\right) = PSD_v \left(f = \frac{1}{2\pi}\right) = PSD_a \left(f = \frac{1}{2\pi}\right), \quad \frac{1}{2\pi} \approx 0.16 [Hz] \]
+
+
+
+
+1.7 Cumulative Power/Amplitude Spectrum
+If we want to compute the Cumulative Power Spectrum: \[ CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2] \] @@ -567,15 +695,425 @@ Then, we can obtain the Root Mean Square value of the velocity:
-
+
++figure; +hold on; +plot(f, cumtrapz(f, Sv)); +hold off; +set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]') ++
+In order to integrate from high frequency to low frequency: +
+
+
+figure; +hold on; +plot(f, flip(-cumtrapz(flip(f), flip(Sv)))); +hold off; +set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]') ++
+
+1.8 TODO Add best practices from the Article and simple snippet that works
+
+
+
+2 Technique 1 : Approximation with a transfer function
+ +
+
+
+2.1 Signal's PSD
+
+
++We load the PSD of the signal we wish to replicate. +
+
+
+
+load('./mat/dist_psd.mat', 'dist_f'); ++
+We remove the first value with very high PSD. +
+
+
+
+dist_f.f = dist_f.f(3:end); +dist_f.psd_gm = dist_f.psd_gm(3:end); ++
+The PSD of the signal is shown on figure fig:psd_ground_motion. +
+ + + +
+
+
+
+2.2 Transfer Function that approximate the ASD
+
+
+
+
+G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2); ++
+
+
+2.3 Generated Time domain signal
+
+
+
+
+
+Fs = 2*dist_f.f(end); +Ts = 1/Fs; + +t = 0:Ts:500; +u = sqrt(Fs/2)*randn(length(t), 1); ++
+
+
+u_o = lsim(G_gm, u, t); ++
+
+figure; +plot(t, u_o); ++
+
+
+2.4 Comparison of the Power Spectral Densities
+
+
+
+
+
+nx = length(u_o); +na = 16; +win = hanning(floor(nx/na)); + +[pxx, f] = pwelch(u_o, win, 0, [], Fs); ++
+Finally, we compare the PSD of the original signal and the obtained signal on figure fig:psd_comparison. +
+
+
2.5 Simulink
+
+
+
+
+
+
+
Figure 8: Simulink Schematic
+
+The parameters for the Band-Limited White Noise are:
+
-
+
- Noise Power: 1 +
+
nx = length(out.u_gm.Data); +na = 8; +win = hanning(floor(nx/na)); + +[pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);
+
+
+3 Technique 2 : IFFT
+
+
+
++The technique comes from preumont94_random_vibrat_spect_analy (section 12.11). +
+
+
+
+3.1 Signal's PSD
+
+
++We load the PSD of the signal we wish to replicate. +
+
+
+
+load('./mat/dist_psd.mat', 'dist_f'); ++
+We remove the first value with very high PSD. +
+
+
+
+dist_f.f = dist_f.f(3:end); +dist_f.psd_gm = dist_f.psd_gm(3:end); ++
+The PSD of the signal is shown on figure fig:psd_original. +
+ + + +
+
+
+
+3.2 Algorithm
+
+
++We define some parameters. +
+
+
+
+Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz] +N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD +T0 = N/Fs; % Signal Duration [s] +df = 1/T0; % Frequency resolution of the DFT [Hz] + % Also equal to (dist_f.f(2)-dist_f.f(1)) ++
+We then specify the wanted PSD. +
+
+
+
+phi = dist_f.psd_gm; ++
+Create amplitudes corresponding to wanted PSD. +
+
+
+
+C = zeros(N/2,1); +for i = 1:N/2 + C(i) = sqrt(phi(i)*df); +end ++
+Add random phase to C.
+
+
+theta = 2*pi*rand(N/2,1); % Generate random phase [rad] + +Cx = [0 ; C.*complex(cos(theta),sin(theta))]; +Cx = [Cx; flipud(conj(Cx(2:end)))];; ++
+
+
+3.3 Obtained Time Domain Signal
+ +
+
+
+3.4 PSD Comparison
+
+
++We duplicate the time domain signal to have a longer signal and thus a more precise PSD result. +
+
+
+
+u_rep = repmat(u, 10, 1); ++
+We compute the PSD of the obtained signal with the following commands. +
+
+
+
+nx = length(u_rep); +na = 16; +win = hanning(floor(nx/na)); + +[pxx, f] = pwelch(u_rep, win, 0, [], Fs); ++
+Finally, we compare the PSD of the original signal and the obtained signal on figure fig:psd_comparison. +
+ + + +
+
+
+
+4 TODO Compute the Noise level and Signal level from PSD
+ +
+
+4.1 Computation
+
+
+-
+
[ ]Add table to compare the methods
+[ ]Add some explanations
+
+
+
+N = 10000; +dt = 0.001; + +t = dt*(0:1:N-1)'; ++
+Parameters of the signal +
+
+
+
+asig = 0.8; % Amplitude of the signal [V] +fsig = 100; % Frequency of the signal [Hz] + +anoi = 1e-3; % RMS value of the noise + +x = anoi*randn(N, 1) + asig*sin((2*pi*fsig)*t); ++
+
+
+figure; +plot(t, x); ++
+Compute the PSD of the signal. +
+
+
+
+nx = length(x); +na = 8; +win = blackmanharris(floor(nx/na)); + +[pxx, f] = pwelch(x, win, 0, [], 1/dt); ++
+Normalization of the PSD. +
+
+
+
+CG = sum(win)/(nx/na); +NG = sum(win.^2)/(nx/na); +fbin = f(2) - f(1); + +pxx_norm = pxx*(NG*fbin/(CG)^2); ++
+
+
+isig = round(fsig/fbin)+1; ++
+Estimate the Signal magnitude. +
+
+
+
+srmt = asig/sqrt(2) % Theoretical value of signal magnitude +srms = sqrt(sum(pxx(isig-5:isig+5)*fbin)) % Signal spectrum integrated +srmsp = sqrt(pxx(isig) * NG*fbin/CG^2) % Maximum read off spectrum ++
+Estimate the noise floor. +
+
+
+
+nth = anoi/sqrt(max(f)) % Theoretical value [V/sqrt(Hz)] + +inmax = isig-20; +nsum = sqrt(sum(pxx(1:inmax)*fbin)) / sqrt(f(inmax)) % Signal spectrum integrated + +navg = sqrt(mean(pxx(1:inmax))) % pwelch output averaged ++
+
+figure; +hold on; +plot(f, pxx) +plot(f, pxx_norm) +hold off; +xlabel('Frequency [Hz]'); +ylabel('Power Spectral Density'); +set(gca, 'xscale', 'log'); +set(gca, 'yscale', 'log'); ++
+ +
Bibliography
+- [preumont94_random_vibrat_spect_analy] Andr\'e Preumont, Random Vibration and Spectral Analysis, Springer Netherlands (1994). +