Update the analysis
BIN
figs/asd_and_tf_compare.png
Normal file
After Width: | Height: | Size: 66 KiB |
BIN
figs/asd_velocity.png
Normal file
After Width: | Height: | Size: 116 KiB |
BIN
figs/compare_psd_original_simulink.png
Normal file
After Width: | Height: | Size: 78 KiB |
BIN
figs/compare_psd_tf_technique.png
Normal file
After Width: | Height: | Size: 78 KiB |
BIN
figs/cps_integral_comp.png
Normal file
After Width: | Height: | Size: 150 KiB |
Before Width: | Height: | Size: 78 KiB After Width: | Height: | Size: 76 KiB |
Before Width: | Height: | Size: 70 KiB After Width: | Height: | Size: 66 KiB |
Before Width: | Height: | Size: 70 KiB After Width: | Height: | Size: 66 KiB |
BIN
figs/psd_velocity.png
Normal file
After Width: | Height: | Size: 144 KiB |
BIN
figs/psd_velocity_displacement.png
Normal file
After Width: | Height: | Size: 82 KiB |
BIN
figs/psd_velocity_lti_method.png
Normal file
After Width: | Height: | Size: 72 KiB |
Before Width: | Height: | Size: 25 KiB After Width: | Height: | Size: 25 KiB |
BIN
figs/time_domain_u.png
Normal file
After Width: | Height: | Size: 32 KiB |
BIN
figs/time_domain_x_zoom.png
Normal file
After Width: | Height: | Size: 64 KiB |
Before Width: | Height: | Size: 47 KiB After Width: | Height: | Size: 37 KiB |
Before Width: | Height: | Size: 7.4 KiB After Width: | Height: | Size: 9.2 KiB |
Before Width: | Height: | Size: 2.4 KiB After Width: | Height: | Size: 2.4 KiB |
Before Width: | Height: | Size: 6.2 KiB After Width: | Height: | Size: 7.8 KiB |
1782
index.html
107
matlab/approximate_psd_ifft.m
Normal file
@ -0,0 +1,107 @@
|
|||||||
|
%% Clear Workspace and Close figures
|
||||||
|
clear; close all; clc;
|
||||||
|
|
||||||
|
%% Intialize Laplace variable
|
||||||
|
s = zpk('s');
|
||||||
|
|
||||||
|
% 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 ref:fig:psd_original.
|
||||||
|
|
||||||
|
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');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
% 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))
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% We then specify the wanted PSD.
|
||||||
|
|
||||||
|
phi = dist_f.psd_gm;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% We create amplitudes corresponding to wanted PSD.
|
||||||
|
|
||||||
|
C = zeros(N/2,1);
|
||||||
|
for i = 1:N/2
|
||||||
|
C(i) = sqrt(phi(i)*df);
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Finally, we add some 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)))];;
|
||||||
|
|
||||||
|
% Obtained Time Domain Signal
|
||||||
|
% 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]
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(t, u)
|
||||||
|
xlabel('Time [s]');
|
||||||
|
ylabel('Amplitude');
|
||||||
|
xlim([t(1), t(end)]);
|
||||||
|
|
||||||
|
% 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 ref:fig:psd_comparison.
|
||||||
|
|
||||||
|
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');
|
||||||
|
xlim([0.1, 500]);
|
150
matlab/approximate_psd_tf.m
Normal file
@ -0,0 +1,150 @@
|
|||||||
|
%% Clear Workspace and Close figures
|
||||||
|
clear; close all; clc;
|
||||||
|
|
||||||
|
%% Intialize Laplace variable
|
||||||
|
s = zpk('s');
|
||||||
|
|
||||||
|
% 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 ref:fig:psd_ground_motion.
|
||||||
|
|
||||||
|
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');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
% 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}$:
|
||||||
|
% \[ |G(j\omega)| \approx \Gamma_x(\omega) \]
|
||||||
|
|
||||||
|
|
||||||
|
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 [[]].
|
||||||
|
|
||||||
|
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('Amplitude Spectral Density');
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
% Generated Time domain signal
|
||||||
|
% We know that a signal $u$ going through a LTI system $G$ (figure [[fig:velocity_to_voltage_psd_bis]]) 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}
|
||||||
|
|
||||||
|
% #+NAME: fig:velocity_to_voltage_psd_bis
|
||||||
|
% #+CAPTION: Schematic of the instrumentation used for the measurement
|
||||||
|
% [[file:figs/velocity_to_voltage_psd.png]]
|
||||||
|
|
||||||
|
% 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 =G_gm=, we should obtain a signal with an ASD that approximate the original ASD.
|
||||||
|
|
||||||
|
% 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]
|
||||||
|
|
||||||
|
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 [[fig:velocity_to_voltage_psd_bis]].
|
||||||
|
|
||||||
|
y = lsim(G_gm, u, t);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% The obtained time domain signal is shown in figure [[fig:time_domain_u]].
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(t, y);
|
||||||
|
xlabel('Time [s]');
|
||||||
|
ylabel('Amplitude');
|
||||||
|
|
||||||
|
% 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));
|
||||||
|
|
||||||
|
[pxx, f] = pwelch(y, win, 0, [], Fs);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% Finally, we compare the PSD of the original signal and the obtained signal on figure ref:fig:psd_comparison.
|
||||||
|
|
||||||
|
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');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
% Simulink
|
||||||
|
% One advantage of this technique is that it can be easily integrated into simulink.
|
||||||
|
|
||||||
|
% The corresponding schematic is shown in figure [[fig:simulink_psd_generate]] 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).
|
||||||
|
|
||||||
|
% Then, the signal generated pass through the transfer function representing the wanted ASD.
|
||||||
|
|
||||||
|
% #+name: fig:simulink_psd_generate
|
||||||
|
% #+caption: Simulink Schematic
|
||||||
|
% [[file:figs/simulink_psd_generate.png]]
|
||||||
|
|
||||||
|
% We simulate the system shown in figure [[fig:simulink_psd_generate]].
|
||||||
|
|
||||||
|
out = sim('matlab/generate_signal_psd.slx');
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% And we compute the PSD of the generated signal.
|
||||||
|
|
||||||
|
nx = length(out.u_gm.Data);
|
||||||
|
na = 8;
|
||||||
|
win = hanning(floor(nx/na));
|
||||||
|
|
||||||
|
[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 [[fig:compare_psd_original_simulink]].
|
||||||
|
|
||||||
|
figure;
|
||||||
|
hold on;
|
||||||
|
plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
|
||||||
|
plot(f, pxx, 'DisplayName', 'Computed from Simulink')
|
||||||
|
hold off;
|
||||||
|
xlabel('Frequency [Hz]');
|
||||||
|
ylabel('Power Spectral Density');
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
legend('location', 'northeast');
|
||||||
|
xlim([0.1, 500]);
|
93
matlab/compute_psd_levels.m
Normal file
@ -0,0 +1,93 @@
|
|||||||
|
%% Clear Workspace and Close figures
|
||||||
|
clear; close all; clc;
|
||||||
|
|
||||||
|
%% Intialize Laplace variable
|
||||||
|
s = zpk('s');
|
||||||
|
|
||||||
|
% 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]
|
||||||
|
|
||||||
|
t = dt*(0:1:N-1)'; % Time vector [s]
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% We generate of signal that consist of:
|
||||||
|
% - a white noise with an RMS value equal to =anoi=
|
||||||
|
% - two sinusoidal signals
|
||||||
|
|
||||||
|
% The parameters are defined below.
|
||||||
|
|
||||||
|
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]
|
||||||
|
|
||||||
|
anoi = 1e-3; % RMS value of the noise
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% The signal $x$ is generated with the following code and is shown in figure [[fig:time_domain_x_zoom]].
|
||||||
|
|
||||||
|
x = anoi*randn(N, 1) + asig*sin((2*pi*fsig)*t) + ahar*sin((2*pi*fhar)*t);
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(t, x);
|
||||||
|
xlabel('Time [s]');
|
||||||
|
ylabel('Amplitude');
|
||||||
|
xlim([0, 1]);
|
||||||
|
|
||||||
|
% 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));
|
||||||
|
|
||||||
|
[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);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% We determine the frequency bins corresponding to the frequency of the signals.
|
||||||
|
|
||||||
|
isig = round(fsig/fbin)+1;
|
||||||
|
ihar = round(fhar/fbin)+1;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% The theoretical RMS value of the signal is:
|
||||||
|
|
||||||
|
srmt = asig/sqrt(2); % Theoretical value of signal magnitude
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% 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
|
||||||
|
|
||||||
|
% Estimation of the noise level
|
||||||
|
% 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)]
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% 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
|
279
matlab/spectral_analysis_basics.m
Normal file
@ -0,0 +1,279 @@
|
|||||||
|
%% Clear Workspace and Close figures
|
||||||
|
clear; close all; clc;
|
||||||
|
|
||||||
|
%% Intialize Laplace variable
|
||||||
|
s = zpk('s');
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% #+RESULTS:
|
||||||
|
% [[file:figs/velocity_to_voltage.png]]
|
||||||
|
|
||||||
|
% #+NAME: fig:velocity_to_voltage
|
||||||
|
% #+CAPTION: Schematic of the instrumentation used for the measurement
|
||||||
|
% #+RESULTS: fig:velocity_to_voltage
|
||||||
|
% [[file:figs/velocity_to_voltage.png]]
|
||||||
|
|
||||||
|
% 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] \]
|
||||||
|
|
||||||
|
% The parameters are defined below.
|
||||||
|
|
||||||
|
G0 = 88; % Sensitivity [V/(m/s)]
|
||||||
|
f0 = 2; % Cut-off frequency [Hz]
|
||||||
|
|
||||||
|
Gg = G0*(s/2/pi/f0)/(1+s/2/pi/f0);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% And the dynamics of the amplifier in the bandwidth of interest is just a gain: $G_a(s) = 1000$.
|
||||||
|
|
||||||
|
Ga = 1000;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% #+NAME: fig:voltage_to_velocity
|
||||||
|
% #+CAPTION: Schematic of the instrumentation used for the measurement
|
||||||
|
% #+RESULTS: fig:voltage_to_velocity
|
||||||
|
% [[file:figs/voltage_to_velocity.png]]
|
||||||
|
|
||||||
|
|
||||||
|
% Let's load the measured $x$.
|
||||||
|
|
||||||
|
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]
|
||||||
|
|
||||||
|
dt = t(2)-t(1); % Sampling Time [s]
|
||||||
|
Fs = 1/dt; % Sampling Frequency [Hz]
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% We simulate this system with matlab using the =lsim= command.
|
||||||
|
|
||||||
|
v = lsim(inv(Gg*Ga), x, t);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% And we plot the obtained velocity
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(t, v);
|
||||||
|
xlabel("Time [s]");
|
||||||
|
ylabel("Velocity [m/s]");
|
||||||
|
|
||||||
|
% Power Spectral Density and Amplitude Spectral Density
|
||||||
|
% From the Matlab documentation:
|
||||||
|
% #+begin_quote
|
||||||
|
% 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.
|
||||||
|
% #+end_quote
|
||||||
|
|
||||||
|
% We now have the velocity $v$ in the time domain:
|
||||||
|
% \[ v(t)\ [m/s] \]
|
||||||
|
|
||||||
|
% 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
|
||||||
|
% - $f_s$ is the sampling frequency in Hz
|
||||||
|
% - $R_{xx}$ is the autocorrelation
|
||||||
|
|
||||||
|
% 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.
|
||||||
|
|
||||||
|
|
||||||
|
% To compute the Power Spectral Density with matlab, we use the =pwelch= function ([[https://fr.mathworks.com/help/signal/ref/pwelch.html?s_tid=doc_ta][documentation]]).
|
||||||
|
% The use of the =pwelch= function is:
|
||||||
|
% =[pxx,w] = pwelch(x,window,noverlap,nfft, fs)=
|
||||||
|
% With:
|
||||||
|
% - =x= is the discrete time signal
|
||||||
|
% - =window= is a window that is used to smooth the obtained PSD
|
||||||
|
% - =overlap= can be used to have some overlap from section to section
|
||||||
|
% - =nfft= specifies the number of FFT points for the PSD
|
||||||
|
% - =fs= is the sampling frequency of the data =x= in Hertz
|
||||||
|
|
||||||
|
% As explained in cite:schmid12_how_to_use_fft_matlab, it is recommended to use the =pwelch= function the following way:
|
||||||
|
|
||||||
|
% First, define a window (preferably the =hanning= one) by specifying the averaging factor =na=.
|
||||||
|
|
||||||
|
nx = length(v);
|
||||||
|
na = 8;
|
||||||
|
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);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% The obtained PSD is shown in figure [[fig:psd_velocity]].
|
||||||
|
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(f, Sv);
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
xlabel('Frequency [Hz]');
|
||||||
|
ylabel('Power Spectral Density $\left[\frac{(m/s)^2}{Hz}\right]$')
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% #+NAME: fig:psd_velocity
|
||||||
|
% #+CAPTION: Power Spectral Density of the measured velocity
|
||||||
|
% #+RESULTS: fig:psd_velocity
|
||||||
|
% [[file:figs/psd_velocity.png]]
|
||||||
|
|
||||||
|
% The Amplitude Spectral Density (ASD) is defined as the square root of the Power Spectral Density and is shown in figure [[fig:asd_velocity]].
|
||||||
|
% \begin{equation}
|
||||||
|
% \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right]
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
|
||||||
|
figure;
|
||||||
|
plot(f, sqrt(Sv));
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
xlabel('Frequency [Hz]');
|
||||||
|
ylabel('Amplitude Spectral Density $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% #+NAME: fig:velocity_to_voltage_psd
|
||||||
|
% #+CAPTION: Schematic of the instrumentation used for the measurement
|
||||||
|
% #+RESULTS: fig:velocity_to_voltage_psd
|
||||||
|
% [[file:figs/velocity_to_voltage_psd.png]]
|
||||||
|
|
||||||
|
% Similarly, the ASD of $y$ is:
|
||||||
|
% \begin{equation}
|
||||||
|
% \Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{uu}(\omega)
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% 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) \]
|
||||||
|
|
||||||
|
% The PSD of $x$ is computed below.
|
||||||
|
|
||||||
|
nx = length(x);
|
||||||
|
na = 8;
|
||||||
|
win = hanning(floor(nx/na));
|
||||||
|
[Sx, f] = pwelch(x, win, 0, [], Fs);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% And the PSD of $v$ is obtained with the below code.
|
||||||
|
|
||||||
|
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 [[fig:psd_velocity_lti_method]].
|
||||||
|
|
||||||
|
|
||||||
|
figure;
|
||||||
|
hold on;
|
||||||
|
plot(f, Sv, 'DisplayName', 'lsim technique');
|
||||||
|
plot(f, Svv, 'DisplayName', 'LTI technique');
|
||||||
|
hold off;
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
xlabel('Frequency [Hz]');
|
||||||
|
ylabel('Power Spectral Density $\left[\frac{(m/s)^2}{Hz}\right]$');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
legend('location', 'southwest');
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% #+NAME: fig:velocity_to_displacement_psd
|
||||||
|
% #+CAPTION: Schematic of the instrumentation used for the measurement
|
||||||
|
% #+RESULTS: fig:velocity_to_displacement_psd
|
||||||
|
% [[file:figs/velocity_to_displacement_psd.png]]
|
||||||
|
|
||||||
|
% We then have the relation between the PSD of $d$ and the PSD of $v$:
|
||||||
|
% \begin{equation}
|
||||||
|
% S_{dd}(\omega) = \left|\frac{1}{j \omega}\right|^2 S_{vv}(\omega)
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% Using a frequency variable $f$ in Hz:
|
||||||
|
% \begin{equation}
|
||||||
|
% S_{dd}(f) = \left| \frac{1}{j 2\pi f} \right|^2 S_{vv}(f)
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% For the Amplitude Spectral Density:
|
||||||
|
% \begin{equation}
|
||||||
|
% \Gamma_{dd}(f) = \frac{1}{2\pi f} \Gamma_{vv}(f)
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% 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$:
|
||||||
|
% \begin{equation}
|
||||||
|
% S_{xx}(\omega = 1) = S_{vv}(\omega = 1)
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
|
||||||
|
% 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;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
% The obtained PSD of the displacement can be seen in figure [[fig:psd_velocity_displacement]].
|
||||||
|
|
||||||
|
figure;
|
||||||
|
hold on;
|
||||||
|
plot(f, Sd, 'DisplayName', '$S_d$ in $\left[\frac{m^2}{Hz}\right]$');
|
||||||
|
plot(f, Sv, 'DisplayName', '$S_v$ in $\left[\frac{(m/s)^2}{Hz}\right]$');
|
||||||
|
hold off;
|
||||||
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||||||
|
xlabel('Frequency [Hz]');
|
||||||
|
ylabel('Power Spectral Density');
|
||||||
|
xlim([0.1, 500]);
|
||||||
|
legend('location', 'southwest');
|
||||||
|
|
||||||
|
% Cumulative Power/Amplitude Spectrum
|
||||||
|
% The Cumulative Power Spectrum is the cumulative integral of the Power Spectral Density:
|
||||||
|
% #+NAME: eq:cps
|
||||||
|
% \begin{equation}
|
||||||
|
% CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2]
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% It is also possible to integrate from high frequency to low frequency:
|
||||||
|
% #+NAME: eq:cps_inv
|
||||||
|
% \begin{equation}
|
||||||
|
% CPS_v(f) = \int_f^\infty PSD_v(\nu) d\nu \quad [(m/s)^2]
|
||||||
|
% \end{equation}
|
||||||
|
|
||||||
|
% 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.
|
||||||
|
|
||||||
|
|
||||||
|
% 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 [[eq:cps_inv]].
|
||||||
|
|
||||||
|
|
||||||
|
% 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] \]
|
||||||
|
|
||||||
|
% 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}] \]
|
||||||
|
|
||||||
|
% With Matlab, the Cumulative Power Spectrum can be computed with the below formulas and the results are shown in figure [[fig:cps_integral_comp]].
|
||||||
|
|
||||||
|
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$]')
|
||||||
|
|
||||||
|
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]);
|
8
ref.bib
@ -18,3 +18,11 @@
|
|||||||
pages = {nil},
|
pages = {nil},
|
||||||
series = {Solid Mechanics and Its Applications},
|
series = {Solid Mechanics and Its Applications},
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@book{stoica05_spect,
|
||||||
|
author = {Stoica, Petre and Moses, Randolph L and others},
|
||||||
|
title = {Spectral analysis of signals},
|
||||||
|
year = 2005,
|
||||||
|
publisher = {Pearson Prentice Hall Upper Saddle River, NJ},
|
||||||
|
keywords = {},
|
||||||
|
}
|