From fffdfd9e91a45b2f47150aea8ea9aa0dbbd599e5 Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Mon, 26 Apr 2021 17:04:09 +0200 Subject: [PATCH] Add link to figure --- index.html | 479 +++++++++++++++++++++++++++-------------------------- index.org | 2 +- 2 files changed, 245 insertions(+), 236 deletions(-) diff --git a/index.html b/index.html index 2427d9c..de9c3b4 100644 --- a/index.html +++ b/index.html @@ -3,21 +3,30 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Spectral Analysis using Matlab - + - - + +
@@ -30,43 +39,43 @@

Table of Contents

@@ -82,23 +91,23 @@ Some matlab documentation about Spectral Analysis can be found 1, some basics of spectral analysis are presented. +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. +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. +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

+
+

1 Spectral Analysis - Basics

- +

In this section, the basics of spectral analysis is presented with the associated Matlab commands. @@ -114,11 +123,11 @@ This include:

-
-

1.1 Sensitivity of the instrumentation

+
+

1.1 Sensitivity of the instrumentation

-A typical measurement setup is shown in figure 1 where we measure a physical signal which is here a velocity \(v(t)\) using a geophone. +A typical measurement setup is shown in figure 1 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)\).

@@ -132,7 +141,7 @@ To obtain the real physical quantity \(v(t)\) as measured by the sensor from the

-
+

velocity_to_voltage.png

Figure 1: Schematic of the instrumentation used for the measurement

@@ -147,10 +156,10 @@ Let’s say, we know that the sensitivity of the geophone used is The parameters are defined below.

-
G0 = 88; % Sensitivity [V/(m/s)]
-f0 = 2; % Cut-off frequency [Hz]
+
  G0 = 88; % Sensitivity [V/(m/s)]
+  f0 = 2; % Cut-off frequency [Hz]
 
-Gg = G0*(s/2/pi/f0)/(1+s/2/pi/f0);
+  Gg = G0*(s/2/pi/f0)/(1+s/2/pi/f0);
@@ -158,21 +167,21 @@ Gg = G0*(s/2
-
Ga = 1000;
+
  Ga = 1000;
-
-

1.2 Convert the time domain from volts to velocity

+
+

1.2 Convert the time domain from volts to velocity

Let’s here try to obtain the time domain signal \(v(t)\) from the measurement \(x\).

-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 2. +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 2.

@@ -180,7 +189,7 @@ If \({G_a(s)}^{-1} {G_g(s)}^{-1}\) is not proper, we add low pass filters at hig

-
+

voltage_to_velocity.png

Figure 2: Schematic of the instrumentation used for the measurement

@@ -191,13 +200,13 @@ If \({G_a(s)}^{-1} {G_g(s)}^{-1}\) is not proper, we add low pass filters at hig Let’s load the measured \(x\).

-
data = load('mat/data_028.mat', 'data'); data = data.data;
+
  data = load('mat/data_028.mat', 'data'); data = data.data;
 
-t = data(:, 3); % Time vector [s]
-x = data(:, 1)-mean(data(:, 1)); % The offset if removed (coming from the voltage amplifier) [v]
+  t = data(:, 3); % Time vector [s]
+  x = data(:, 1)-mean(data(:, 1)); % The offset if removed (coming from the voltage amplifier) [v]
 
-dt = t(2)-t(1); % Sampling Time [s]
-Fs = 1/dt; % Sampling Frequency [Hz]
+  dt = t(2)-t(1); % Sampling Time [s]
+  Fs = 1/dt; % Sampling Frequency [Hz]
@@ -205,7 +214,7 @@ Fs = 1/dt; % Sampling Fr We simulate this system with matlab using the lsim command.

-
v = lsim(inv(Gg*Ga), x, t);
+
  v = lsim(inv(Gg*Ga), x, t);
@@ -213,7 +222,7 @@ We simulate this system with matlab using the lsim command. And we plot the obtained velocity

-
+

velocity_time.png

Figure 3: Computed Velocity from the measured Voltage

@@ -221,8 +230,8 @@ And we plot the obtained velocity
-
-

1.3 Power Spectral Density and Amplitude Spectral Density

+
+

1.3 Power Spectral Density and Amplitude Spectral Density

From the Matlab documentation: @@ -275,9 +284,9 @@ As explained in (Schmid 2012), it is recommen First, define a window (preferably the hanning one) by specifying the averaging factor na.

-
nx = length(v);
-na = 8;
-win = hanning(floor(nx/na));
+
  nx = length(v);
+  na = 8;
+  win = hanning(floor(nx/na));
@@ -285,30 +294,30 @@ win = hanning(floor(nx/na)); Then, compute the power spectral density \(S_v\) and the associated frequency vector \(f\) with no overlap.

-
[Sv, f] = pwelch(v, win, 0, [], Fs);
+
  [Sv, f] = pwelch(v, win, 0, [], Fs);

-The obtained PSD is shown in figure 4. +The obtained PSD is shown in figure 4.

-
+

psd_velocity.png

Figure 4: Power Spectral Density of the measured velocity

-The Amplitude Spectral Density (ASD) is defined as the square root of the Power Spectral Density and is shown in figure 5. +The Amplitude Spectral Density (ASD) is defined as the square root of the Power Spectral Density and is shown in figure 5.

\begin{equation} \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right] \end{equation} -
+

asd_velocity.png

Figure 5: Power Spectral Density of the measured velocity

@@ -316,22 +325,22 @@ The Amplitude Spectral Density (ASD) is defined as the square root of the Power
-
-

1.4 Modification of a signal’s Power Spectral Density when going through an LTI system

+
+

1.4 Modification of a signal’s Power Spectral Density when going through an LTI system

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.

-To do so, we use the fact that a signal \(u\) with a PSD \(S_{uu}\) going through a LTI system \(G_(s)\) (figure 6) will generate a signal \(y\) with a PSD: +To do so, we use the fact that a signal \(u\) with a PSD \(S_{uu}\) going through a LTI system \(G_(s)\) (figure 6) will generate a signal \(y\) with a PSD:

\begin{equation} S_{yy}(\omega) = \left|G(j\omega)\right|^2 S_{uu}(\omega) \end{equation} -
+

velocity_to_voltage_psd.png

Figure 6: Schematic of the instrumentation used for the measurement

@@ -353,10 +362,10 @@ Thus, we could have computed the PSD of \(x\) and then obtain the PSD of the vel The PSD of \(x\) is computed below.

-
nx = length(x);
-na = 8;
-win = hanning(floor(nx/na));
-[Sx, f] = pwelch(x, win, 0, [], Fs);
+
  nx = length(x);
+  na = 8;
+  win = hanning(floor(nx/na));
+  [Sx, f] = pwelch(x, win, 0, [], Fs);
@@ -364,16 +373,16 @@ win = hanning(floor(nx/na)); And the PSD of \(v\) is obtained with the below code.

-
Svv = Sx.*squeeze(abs(freqresp(inv(Gg*Ga), f, 'Hz'))).^2;
+
  Svv = Sx.*squeeze(abs(freqresp(inv(Gg*Ga), f, 'Hz'))).^2;

-The result is compare with the PSD computed from the \(v\) signal obtained with the lsim command in figure 7. +The result is compare with the PSD computed from the \(v\) signal obtained with the lsim command in figure 7.

-
+

psd_velocity_lti_method.png

Figure 7: Obtain PSD of the velocity using the formula (png, pdf)

@@ -381,15 +390,15 @@ The result is compare with the PSD computed from the \(v\) signal obtained with
-
-

1.5 From PSD of the velocity to the PSD of the displacement

+
+

1.5 From PSD of the velocity to the PSD of the displacement

-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 8. +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 8.

-
+

velocity_to_displacement_psd.png

Figure 8: Schematic of the instrumentation used for the measurement

@@ -428,15 +437,15 @@ Note here that the PSD (and ASD) of one variable and its derivatives/integrals a With Matlab, the PSD of the displacement can be computed from the PSD of the velocity with the following code.

-
Sd = Sv.*(1./(2*pi*f)).^2;
+
  Sd = Sv.*(1./(2*pi*f)).^2;

-The obtained PSD of the displacement can be seen in figure 9. +The obtained PSD of the displacement can be seen in figure 9.

-
+

psd_velocity_displacement.png

Figure 9: PSD of the Velocity and Displacement (png, pdf)

@@ -444,14 +453,14 @@ The obtained PSD of the displacement can be seen in figure
-
-

1.6 Cumulative Power/Amplitude Spectrum

+
+

1.6 Cumulative Power/Amplitude Spectrum

The Cumulative Power Spectrum is the cumulative integral of the Power Spectral Density:

\begin{equation} -\label{orgf27b554} +\label{org644886e} CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2] \end{equation} @@ -459,7 +468,7 @@ The Cumulative Power Spectrum is the cumulative integral of the Power Spectral D It is also possible to integrate from high frequency to low frequency:

\begin{equation} -\label{org1c25e5d} +\label{org47d268d} CPS_v(f) = \int_f^\infty PSD_v(\nu) d\nu \quad [(m/s)^2] \end{equation} @@ -470,7 +479,7 @@ The Cumulative Power Spectrum taken at frequency \(f\) thus represent the power

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{org1c25e5d}. +If the power is mostly present at low frequencies, it is preferable to use equation \eqref{org47d268d}.

@@ -485,37 +494,37 @@ The Root Mean Square value of the velocity corresponds to the Cumulative Amplitu

-With Matlab, the Cumulative Power Spectrum can be computed with the below formulas and the results are shown in figure 10. +With Matlab, the Cumulative Power Spectrum can be computed with the below formulas and the results are shown in figure 10.

-
CPS_v  = cumtrapz(f, Sv); % Cumulative Power Spectrum from low to high frequencies
-CPS_vv = flip(-cumtrapz(flip(f), flip(Sv))); % Cumulative Power Spectrum from high to low frequencies
+
  CPS_v  = cumtrapz(f, Sv); % Cumulative Power Spectrum from low to high frequencies
+  CPS_vv = flip(-cumtrapz(flip(f), flip(Sv))); % Cumulative Power Spectrum from high to low frequencies
-
figure;
-ax1 = subplot(1, 2, 1);
-hold on;
-plot(f, CPS_v);
-hold off;
-set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
-xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]')
+
  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(f, CPS_v);
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum [$(m/s)^2$]')
 
-ax2 = subplot(1, 2, 2);
-hold on;
-plot(f, CPS_vv);
-hold off;
-set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
-xlabel('Frequency [Hz]');
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(f, CPS_vv);
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]');
 
-linkaxes([ax1,ax2],'xy');
-xlim([0.1, 500]); ylim([1e-15, 1e-11]);
+  linkaxes([ax1,ax2],'xy');
+  xlim([0.1, 500]); ylim([1e-15, 1e-11]);
-
+

cps_integral_comp.png

Figure 10: Cumulative Power Spectrum (png, pdf)

@@ -524,22 +533,22 @@ xlim([0.1, 500]); ylim([1e-15, 1e -

2 Time domain signal that approximate a PSD - TF technique

+
+

2 Time domain signal that approximate a PSD - TF technique

- +

-
-

2.1 Signal’s PSD

+
+

2.1 Signal’s PSD

We load the PSD of the signal we wish to replicate.

-
load('./mat/dist_psd.mat', 'dist_f');
+
  load('./mat/dist_psd.mat', 'dist_f');
@@ -547,8 +556,8 @@ We load the PSD of the signal we wish to replicate. 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);
+
  dist_f.f = dist_f.f(3:end);
+  dist_f.psd_gm = dist_f.psd_gm(3:end);
@@ -556,7 +565,7 @@ dist_f.psd_gm = dist_f.psd_gm(3:end); The PSD of the signal is shown on figure fig:psd_ground_motion.

-
+

psd_ground_motion.png

Figure 11: PSD of the signal (png, pdf)

@@ -564,8 +573,8 @@ The PSD of the signal is shown on figure fig:ps
-
-

2.2 Transfer Function that approximate the ASD

+
+

2.2 Transfer Function that approximate the ASD

Using sisotool 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}\): @@ -573,15 +582,15 @@ Using sisotool or any other tool, we create a transfer function \(G

-
G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2);
+
  G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2);

-We compare the ASD \(\Gamma_x(\omega)\) and the magnitude of the generated transfer function \(|G(j\omega)|\) in figure [[]]. +We compare the ASD \(\Gamma_x(\omega)\) and the magnitude of the generated transfer function \(|G(j\omega)|\) in figure 12.

-
+

asd_and_tf_compare.png

Figure 12: Comparison of the ASD and of the transfer function’s magnitude (png, pdf)

@@ -589,18 +598,18 @@ We compare the ASD \(\Gamma_x(\omega)\) and the magnitude of the generated trans
-
-

2.3 Generated Time domain signal

+
+

2.3 Generated Time domain signal

-We know that a signal \(u\) going through a LTI system \(G\) (figure 13) will have its ASD modified according to the following equation: +We know that a signal \(u\) going through a LTI system \(G\) (figure 13) will have its ASD modified according to the following equation:

\begin{equation} \Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{uu}(\omega) \end{equation} -
+

velocity_to_voltage_psd.png

Figure 13: Schematic of the instrumentation used for the measurement

@@ -614,27 +623,27 @@ Thus, if we create a random signal with an ASD equal to one at all frequency and To obtain a random signal with an ASD equal to one, we use the following code.

-
Fs = 2*dist_f.f(end); % Sampling Frequency [Hz]
-Ts = 1/Fs; % Sampling Time [s]
+
  Fs = 2*dist_f.f(end); % Sampling Frequency [Hz]
+  Ts = 1/Fs; % Sampling Time [s]
 
-t = 0:Ts:500; % Time Vector [s]
-u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one
+  t = 0:Ts:500; % Time Vector [s]
+  u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one

-We then use lsim to compute \(y\) as shown in figure 13. +We then use lsim to compute \(y\) as shown in figure 13.

-
y = lsim(G_gm, u, t);
+
  y = lsim(G_gm, u, t);

-The obtained time domain signal is shown in figure 14. +The obtained time domain signal is shown in figure 14.

-
+

time_domain_u.png

Figure 14: Obtained time domain signal \(y(t)\) (png, pdf)

@@ -642,18 +651,18 @@ The obtained time domain signal is shown in figure 14.
-
-

2.4 Comparison of the Power Spectral Densities

+
+

2.4 Comparison of the Power Spectral Densities

We now compute the Power Spectral Density of the computed time domain signal \(y\).

-
nx = length(y);
-na = 16;
-win = hanning(floor(nx/na));
+
  nx = length(y);
+  na = 16;
+  win = hanning(floor(nx/na));
 
-[pxx, f] = pwelch(y, win, 0, [], Fs);
+  [pxx, f] = pwelch(y, win, 0, [], Fs);
@@ -661,7 +670,7 @@ win = hanning(floor(nx/na)); Finally, we compare the PSD of the original signal and the obtained signal on figure fig:psd_comparison.

-
+

compare_psd_tf_technique.png

Figure 15: Comparison of the original PSD and the PSD of the computed time domain signal (png, pdf)

@@ -669,15 +678,15 @@ Finally, we compare the PSD of the original signal and the obtained signal on fi
-
-

2.5 Simulink

+
+

2.5 Simulink

One advantage of this technique is that it can be easily integrated into simulink.

-The corresponding schematic is shown in figure 16 where the block Band-Limited White Noise is used to generate a random signal with a PSD equal to one (parameter Noise Power is set to 1). +The corresponding schematic is shown in figure 16 where the block Band-Limited White Noise is used to generate a random signal with a PSD equal to one (parameter Noise Power is set to 1).

@@ -685,17 +694,17 @@ Then, the signal generated pass through the transfer function representing the w

-
+

simulink_psd_generate.png

Figure 16: Simulink Schematic

-We simulate the system shown in figure 16. +We simulate the system shown in figure 16.

-
out = sim('matlab/generate_signal_psd.slx');
+
  out = sim('matlab/generate_signal_psd.slx');
@@ -703,19 +712,19 @@ We simulate the system shown in figure 16. And we compute the PSD of the generated signal.

-
nx = length(out.u_gm.Data);
-na = 8;
-win = hanning(floor(nx/na));
+
  nx = length(out.u_gm.Data);
+  na = 8;
+  win = hanning(floor(nx/na));
 
-[pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);
+  [pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);

-Finally, we compare the PSD of the generated signal with the original PSD in figure 17. +Finally, we compare the PSD of the generated signal with the original PSD in figure 17.

-
+

compare_psd_original_simulink.png

Figure 17: Comparison of the obtained signal’s PSD and original PSD (png, pdf)

@@ -724,11 +733,11 @@ Finally, we compare the PSD of the generated signal with the original PSD in fig
-
-

3 Time domain signal that approximate a PSD - IFFT technique

+
+

3 Time domain signal that approximate a PSD - IFFT technique

- +

The technique comes from (Preumont 1994) (section 12.11). @@ -737,14 +746,14 @@ It makes used of the Unversed Fast Fourier Transform (IFFT).

-
-

3.1 Signal’s PSD

+
+

3.1 Signal’s PSD

We load the PSD of the signal we wish to replicate.

-
load('./mat/dist_psd.mat', 'dist_f');
+
  load('./mat/dist_psd.mat', 'dist_f');
@@ -752,8 +761,8 @@ We load the PSD of the signal we wish to replicate. 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);
+
  dist_f.f = dist_f.f(3:end);
+  dist_f.psd_gm = dist_f.psd_gm(3:end);
@@ -761,7 +770,7 @@ dist_f.psd_gm = dist_f.psd_gm(3:end); The PSD of the signal is shown on figure fig:psd_original.

-
+

psd_original.png

Figure 18: PSD of the original signal (png, pdf)

@@ -769,18 +778,18 @@ The PSD of the signal is shown on figure fig:psd_ori
-
-

3.2 Algorithm

+
+

3.2 Algorithm

We define some parameters that will be used in the algorithm.

-
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))
+
  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))
@@ -788,7 +797,7 @@ df = 1/T0;
-
phi = dist_f.psd_gm;
+
  phi = dist_f.psd_gm;
@@ -796,10 +805,10 @@ We then specify the wanted PSD. We create amplitudes corresponding to wanted PSD.

-
C = zeros(N/2,1);
-for i = 1:N/2
-  C(i) = sqrt(phi(i)*df);
-end
+
  C = zeros(N/2,1);
+  for i = 1:N/2
+    C(i) = sqrt(phi(i)*df);
+  end
@@ -807,17 +816,17 @@ We create amplitudes corresponding to wanted PSD. Finally, we add some random phase to C.

-
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+
  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)))];;
+  Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+  Cx = [Cx; flipud(conj(Cx(2:end)))];;
-
-

3.3 Obtained Time Domain Signal

+
+

3.3 Obtained Time Domain Signal

The time domain data is generated by an inverse FFT. @@ -827,13 +836,13 @@ The time domain data is generated by an inverse FFT. The ifft Matlab does not take into account the sampling frequency, thus we need to normalize the signal.

-
u = N/sqrt(2)*ifft(Cx); % Normalisation of the IFFT
-t = linspace(0, T0, N+1); % Time Vector [s]
+
  u = N/sqrt(2)*ifft(Cx); % Normalisation of the IFFT
+  t = linspace(0, T0, N+1); % Time Vector [s]
-
+

signal_time_domain.png

Figure 19: Obtained signal in the time domain (png, pdf)

@@ -841,14 +850,14 @@ t = linspace(0, T0, N+1); -

3.4 PSD Comparison

+
+

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);
+
  u_rep = repmat(u, 10, 1);
@@ -856,11 +865,11 @@ We duplicate the time domain signal to have a longer signal and thus a more prec We compute the PSD of the obtained signal with the following commands.

-
nx = length(u_rep);
-na = 16;
-win = hanning(floor(nx/na));
+
  nx = length(u_rep);
+  na = 16;
+  win = hanning(floor(nx/na));
 
-[pxx, f] = pwelch(u_rep, win, 0, [], Fs);
+  [pxx, f] = pwelch(u_rep, win, 0, [], Fs);
@@ -868,7 +877,7 @@ win = hanning(floor(nx/na)); Finally, we compare the PSD of the original signal and the obtained signal on figure fig:psd_comparison.

-
+

psd_comparison.png

Figure 20: Comparison of the PSD of the original signal and the PSD of the obtained signal (png, pdf)

@@ -877,11 +886,11 @@ Finally, we compare the PSD of the original signal and the obtained signal on fi
-
-

4 Compute the Noise level and Signal level from PSD

+
+

4 Compute the Noise level and Signal level from PSD

- +

We here make use of the Power Spectral Density to estimate either the noise level or the amplitude of a deterministic signal. @@ -889,17 +898,17 @@ Everything is explained in (Schmid 2012) sect

-
-

4.1 Time Domain Signal

+
+

4.1 Time Domain Signal

Let’s first define the number of sample and the sampling time.

-
N  = 10000; % Number of Sample
-dt = 0.001; % Sampling Time [s]
+
  N  = 10000; % Number of Sample
+  dt = 0.001; % Sampling Time [s]
 
-t = dt*(0:1:N-1)'; % Time vector [s]
+  t = dt*(0:1:N-1)'; % Time vector [s]
@@ -915,26 +924,26 @@ We generate of signal that consist of: The parameters are defined below.

-
asig = 0.1; % Amplitude of the signal [V]
-fsig = 10;  % Frequency of the signal [Hz]
+
  asig = 0.1; % Amplitude of the signal [V]
+  fsig = 10;  % Frequency of the signal [Hz]
 
-ahar = 0.5; % Amplitude of the harmonic [V]
-fhar = 50;  % Frequency of the harmonic [Hz]
+  ahar = 0.5; % Amplitude of the harmonic [V]
+  fhar = 50;  % Frequency of the harmonic [Hz]
 
-anoi = 1e-3; % RMS value of the noise
+  anoi = 1e-3; % RMS value of the noise

-The signal \(x\) is generated with the following code and is shown in figure 21. +The signal \(x\) is generated with the following code and is shown in figure 21.

-
x = anoi*randn(N, 1) + asig*sin((2*pi*fsig)*t) + ahar*sin((2*pi*fhar)*t);
+
  x = anoi*randn(N, 1) + asig*sin((2*pi*fsig)*t) + ahar*sin((2*pi*fhar)*t);
-
+

time_domain_x_zoom.png

Figure 21: Time Domain Signal (png, pdf)

@@ -942,18 +951,18 @@ The signal \(x\) is generated with the following code and is shown in figure
-
-

4.2 Estimation of the magnitude of a deterministic signal

+
+

4.2 Estimation of the magnitude of a deterministic signal

Let’s compute the PSD of the signal using the blackmanharris window.

-
nx = length(x);
-na = 8;
-win = blackmanharris(floor(nx/na));
+
  nx = length(x);
+  na = 8;
+  win = blackmanharris(floor(nx/na));
 
-[pxx, f] = pwelch(x, win, 0, [], 1/dt);
+  [pxx, f] = pwelch(x, win, 0, [], 1/dt);
@@ -961,11 +970,11 @@ win = blackmanharris(floor(nx/na)); Normalization of the PSD.

-
CG = sum(win)/(nx/na);
-NG = sum(win.^2)/(nx/na);
-fbin = f(2) - f(1);
+
  CG = sum(win)/(nx/na);
+  NG = sum(win.^2)/(nx/na);
+  fbin = f(2) - f(1);
 
-pxx_norm = pxx*(NG*fbin/CG^2);
+  pxx_norm = pxx*(NG*fbin/CG^2);
@@ -973,8 +982,8 @@ pxx_norm = pxx*(NG*f We determine the frequency bins corresponding to the frequency of the signals.

-
isig = round(fsig/fbin)+1;
-ihar = round(fhar/fbin)+1;
+
  isig = round(fsig/fbin)+1;
+  ihar = round(fhar/fbin)+1;
@@ -982,7 +991,7 @@ ihar = round(fhar/fbin)+
-
srmt = asig/sqrt(2); % Theoretical value of signal magnitude
+
  srmt = asig/sqrt(2); % Theoretical value of signal magnitude
@@ -990,8 +999,8 @@ The theoretical RMS value of the signal is: 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.

-
srms = sqrt(sum(pxx(isig-5:isig+5)*fbin)); % Signal spectrum integrated
-srmsp = sqrt(pxx_norm(isig) * NG*fbin/CG^2); % Maximum read off spectrum
+
  srms = sqrt(sum(pxx(isig-5:isig+5)*fbin)); % Signal spectrum integrated
+  srmsp = sqrt(pxx_norm(isig) * NG*fbin/CG^2); % Maximum read off spectrum
@@ -1035,8 +1044,8 @@ Thus, always the integrated method should be used.
-
-

4.3 Estimation of the noise level

+
+

4.3 Estimation of the noise level

The noise level can also be computed using the integration method. @@ -1046,7 +1055,7 @@ The noise level can also be computed using the integration method. The theoretical RMS noise value is.

-
nth = anoi/sqrt(max(f)) % Theoretical value [V/sqrt(Hz)]
+
  nth = anoi/sqrt(max(f)) % Theoretical value [V/sqrt(Hz)]
@@ -1054,7 +1063,7 @@ The theoretical RMS noise value is. 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.

-
navg = sqrt(mean(pxx_norm([ihar+10:end]))) % pwelch output averaged
+
  navg = sqrt(mean(pxx_norm([ihar+10:end]))) % pwelch output averaged
@@ -1092,12 +1101,12 @@ The estimate of the noise level is quite good.
-
-

5 Further Notes

+
+

5 Further Notes

-
-

5.1 PSD of ADC quantization noise

+
+

5.1 PSD of ADC quantization noise

This is taken from here. @@ -1121,11 +1130,11 @@ 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. 22). +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. 22).

-
+

probability_density_function_adc.png

Figure 22: Probability density function \(p(e)\) of the ADC error \(e\)

@@ -1183,7 +1192,7 @@ Finally:

Author: Dehaeze Thomas

-

Created: 2020-11-12 jeu. 10:22

+

Created: 2021-04-26 lun. 17:04

diff --git a/index.org b/index.org index 0fb5a1b..a107478 100644 --- a/index.org +++ b/index.org @@ -533,7 +533,7 @@ Using =sisotool= or any other tool, we create a transfer function $G$ such that G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2); #+end_src -We compare the ASD $\Gamma_x(\omega)$ and the magnitude of the generated transfer function $|G(j\omega)|$ in figure [[]]. +We compare the ASD $\Gamma_x(\omega)$ and the magnitude of the generated transfer function $|G(j\omega)|$ in figure [[fig:asd_and_tf_compare]]. #+begin_src matlab :exports none figure; hold on;