153 lines
4.3 KiB
Matlab
153 lines
4.3 KiB
Matlab
%% Clear Workspace and Close figures
|
|
clear; close all; clc;
|
|
|
|
%% Intialize Laplace variable
|
|
s = zpk('s');
|
|
|
|
addpath('./mat/');
|
|
|
|
% 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]);
|