Finish the file about measurements. Add raw measurements to git
This commit is contained in:
317
modal-analysis/matlab/modal_frf_coh.m
Normal file
317
modal-analysis/matlab/modal_frf_coh.m
Normal file
@@ -0,0 +1,317 @@
|
||||
%% Clear Workspace and Close figures
|
||||
clear; close all; clc;
|
||||
|
||||
%% Intialize Laplace variable
|
||||
s = zpk('s');
|
||||
|
||||
% Windowing
|
||||
% Windowing is used on the force and response signals.
|
||||
|
||||
% A boxcar window (figure [[fig:windowing_force_signal]]) is used for the force signal as once the impact on the structure is done, the measured signal is meaningless.
|
||||
% The parameters are:
|
||||
% - *Start*: $3\%$
|
||||
% - *Stop*: $7\%$
|
||||
|
||||
|
||||
figure;
|
||||
plot(100*[0, 0.03, 0.03, 0.07, 0.07, 1], [0, 0, 1, 1, 0, 0]);
|
||||
xlabel('Time [%]'); ylabel('Amplitude');
|
||||
xlim([0, 100]); ylim([0, 1]);
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:windowing_force_signal
|
||||
% #+CAPTION: Window used for the force signal
|
||||
% [[file:figs/windowing_force_signal.png]]
|
||||
|
||||
% An exponential window (figure [[fig:windowing_response_signal]]) is used for the response signal as we are measuring transient signals and most of the information is located at the beginning of the signal.
|
||||
% The parameters are:
|
||||
% - FlatTop:
|
||||
% - *Start*: $3\%$
|
||||
% - *Stop*: $2.96\%$
|
||||
% - Decreasing point:
|
||||
% - *X*: $60.4\%$
|
||||
% - *Y*: $14.7\%$
|
||||
|
||||
|
||||
x0 = 0.296;
|
||||
xd = 0.604;
|
||||
yd = 0.147;
|
||||
|
||||
alpha = log(yd)/(x0 - xd);
|
||||
|
||||
t = x0:0.01:1.01;
|
||||
y = exp(-alpha*(t-x0));
|
||||
|
||||
figure;
|
||||
plot(100*[0, 0.03, 0.03, x0, t], [0, 0, 1, 1, y]);
|
||||
xlabel('Time [%]'); ylabel('Amplitude');
|
||||
xlim([0, 100]); ylim([0, 1]);
|
||||
|
||||
% Force and Response signals
|
||||
% Let's load some obtained data to look at the force and response signals.
|
||||
|
||||
|
||||
meas1_raw = load('mat/meas_raw_1.mat');
|
||||
|
||||
% Raw Force Data
|
||||
% The force input for the first measurement is shown on figure [[fig:raw_data_force]]. We can see 10 impacts, one zoom on one impact is shown on figure [[fig:raw_data_force_zoom]].
|
||||
|
||||
% The Fourier transform of the force is shown on figure [[fig:fourier_transfor_force_impact]]. This has been obtained without any windowing.
|
||||
|
||||
|
||||
time = linspace(0, meas1_raw.Track1_X_Resolution*length(meas1_raw.Track1), length(meas1_raw.Track1));
|
||||
|
||||
figure;
|
||||
plot(time, meas1_raw.Track1);
|
||||
xlabel('Time [s]');
|
||||
ylabel('Force [N]');
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:raw_data_force
|
||||
% #+CAPTION: Raw Force Data from Hammer Blow
|
||||
% [[file:figs/raw_data_force.png]]
|
||||
|
||||
|
||||
figure;
|
||||
plot(time, meas1_raw.Track1);
|
||||
xlabel('Time [s]');
|
||||
ylabel('Force [N]');
|
||||
xlim([22.1, 22.3]);
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:raw_data_force_zoom
|
||||
% #+CAPTION: Raw Force Data from Hammer Blow - Zoom
|
||||
% [[file:figs/raw_data_force_zoom.png]]
|
||||
|
||||
|
||||
Fs = 1/meas1_raw.Track1_X_Resolution; % Sampling Frequency [Hz]
|
||||
impacts = [5.9, 11.2, 16.6, 22.2, 27.3, 32.7, 38.1, 43.8, 50.4]; % Time just before the impact occurs [s]
|
||||
|
||||
L = 8194;
|
||||
f = Fs*(0:(L/2))/L; % Frequency vector [Hz]
|
||||
|
||||
F_fft = zeros((8193+1)/2+1, length(impacts));
|
||||
|
||||
for i = 1:length(impacts)
|
||||
t0 = impacts(i);
|
||||
[~, i_start] = min(abs(time-t0));
|
||||
i_end = i_start + 8193;
|
||||
|
||||
Y = fft(meas1_raw.Track1(i_start:i_end));
|
||||
|
||||
P2 = abs(Y/L);
|
||||
P1 = P2(1:L/2+1);
|
||||
P1(2:end-1) = 2*P1(2:end-1);
|
||||
F_fft(:, i) = P1;
|
||||
end
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
for i = 1:length(impacts)
|
||||
plot(f, F_fft(:, i), '-k');
|
||||
end
|
||||
hold off;
|
||||
xlim([0, 200]);
|
||||
xlabel('Frequency [Hz]'); ylabel('Force [N]');
|
||||
|
||||
% Raw Response Data
|
||||
% The response signal for the first measurement is shown on figure [[fig:raw_data_acceleration]]. One zoom on one response is shown on figure [[fig:raw_data_acceleration_zoom]].
|
||||
|
||||
% The Fourier transform of the response signals is shown on figure [[fig:fourier_transform_response_signals]]. This has been obtained without any windowing.
|
||||
|
||||
|
||||
figure;
|
||||
plot(time, meas1_raw.Track2);
|
||||
xlabel('Time [s]');
|
||||
ylabel('Acceleration [m/s2]');
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:raw_data_acceleration
|
||||
% #+CAPTION: Raw Acceleration Data from Accelerometer
|
||||
% [[file:figs/raw_data_acceleration.png]]
|
||||
|
||||
|
||||
figure;
|
||||
plot(time, meas1_raw.Track2);
|
||||
xlabel('Time [s]');
|
||||
ylabel('Acceleration [m/s2]');
|
||||
xlim([22.1, 22.5]);
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:raw_data_acceleration_zoom
|
||||
% #+CAPTION: Raw Acceleration Data from Accelerometer - Zoom
|
||||
% [[file:figs/raw_data_acceleration_zoom.png]]
|
||||
|
||||
|
||||
X_fft = zeros((8193+1)/2+1, length(impacts));
|
||||
|
||||
for i = 1:length(impacts)
|
||||
t0 = impacts(i);
|
||||
[~, i_start] = min(abs(time-t0));
|
||||
i_end = i_start + 8193;
|
||||
|
||||
Y = fft(meas1_raw.Track2(i_start:i_end));
|
||||
|
||||
P2 = abs(Y/L);
|
||||
P1 = P2(1:L/2+1);
|
||||
P1(2:end-1) = 2*P1(2:end-1);
|
||||
X_fft(:, i) = P1;
|
||||
end
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
for i = 1:length(impacts)
|
||||
plot(f, X_fft(:, i), '-k');
|
||||
end
|
||||
hold off;
|
||||
xlim([0, 200]);
|
||||
set(gca, 'Yscale', 'log');
|
||||
xlabel('Frequency [Hz]'); ylabel('Acceleration [$m/s^2$]');
|
||||
|
||||
% Computation of one Frequency Response Function
|
||||
% Now that we have obtained the Fourier transform of both the force input and the response signal, we can compute the Frequency Response Function from the force to the acceleration.
|
||||
|
||||
% We then compare the result obtained with the FRF computed by the modal software (figure [[fig:frf_comparison_software]]).
|
||||
|
||||
% The slight difference can probably be explained by the use of windows.
|
||||
|
||||
% In the following analysis, FRF computed from the software will be used.
|
||||
|
||||
|
||||
meas1 = load('mat/meas_frf_coh_1.mat');
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
for i = 1:length(impacts)
|
||||
plot(f, X_fft(:, i)./F_fft(:, i), '-k');
|
||||
end
|
||||
plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_AvSpc_2_RMS_Y_Val)
|
||||
hold off;
|
||||
xlim([0, 200]);
|
||||
set(gca, 'Yscale', 'log');
|
||||
xlabel('Frequency [Hz]'); ylabel('Acceleration/Force [$m/s^2/N$]');
|
||||
|
||||
% Frequency Response Functions and Coherence Results
|
||||
% Let's see one computed Frequency Response Function and one coherence in order to attest the quality of the measurement.
|
||||
|
||||
% First, we load the data.
|
||||
|
||||
meas1 = load('mat/meas_frf_coh_1.mat');
|
||||
|
||||
|
||||
|
||||
% And we plot on figure [[fig:frf_result_example]] the frequency response function from the force applied in the $X$ direction at the location of the accelerometer number 11 to the acceleration in the $X$ direction as measured by the first accelerometer located on the top platform of the hexapod.
|
||||
|
||||
% The coherence associated is shown on figure [[fig:frf_result_example]].
|
||||
|
||||
|
||||
figure;
|
||||
ax1 = subplot(2, 1, 1);
|
||||
plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_AvXSpc_2_1_RMS_Y_Mod);
|
||||
set(gca, 'Yscale', 'log');
|
||||
set(gca, 'XTickLabel',[]);
|
||||
ylabel('Magnitude');
|
||||
|
||||
ax2 = subplot(2, 1, 2);
|
||||
plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_AvXSpc_2_1_RMS_Y_Phas);
|
||||
ylim([-180, 180]);
|
||||
yticks([-180, -90, 0, 90, 180]);
|
||||
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
||||
|
||||
linkaxes([ax1,ax2],'x');
|
||||
|
||||
|
||||
|
||||
% #+NAME: fig:frf_result_example
|
||||
% #+CAPTION: Example of one measured FRF
|
||||
% [[file:figs/frf_result_example.png]]
|
||||
|
||||
|
||||
figure;
|
||||
plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_Coh_2_1_RMS_Y_Val);
|
||||
xlabel('Frequency [Hz]');
|
||||
ylabel('Coherence');
|
||||
|
||||
% Generation of a FRF matrix and a Coherence matrix from the measurements
|
||||
% We want here to combine all the Frequency Response Functions measured into one big array called the *Frequency Response Matrix*.
|
||||
|
||||
% The frequency response matrix is an $n \times p \times q$:
|
||||
% - $n$ is the number of measurements: $23 \times 3$ (23 accelerometers measuring 3 directions each)
|
||||
% - $p$ is the number of excitation inputs: $3$
|
||||
% - $q$ is the number of frequency points $\omega_i$
|
||||
|
||||
% Thus, the FRF matrix is an $69 \times 3 \times 801$ matrix.
|
||||
|
||||
|
||||
% #+begin_important
|
||||
% For each frequency point $\omega_i$, we obtain a 2D matrix:
|
||||
% \begin{equation}
|
||||
% \text{FRF}(\omega_i) = \begin{bmatrix}
|
||||
% \frac{D_{1_x}}{F_x}(\omega_i) & \frac{D_{1_x}}{F_y}(\omega_i) & \frac{D_{1_x}}{F_z}(\omega_i) \\
|
||||
% \frac{D_{1_y}}{F_x}(\omega_i) & \frac{D_{1_y}}{F_y}(\omega_i) & \frac{D_{1_y}}{F_z}(\omega_i) \\
|
||||
% \frac{D_{1_z}}{F_x}(\omega_i) & \frac{D_{1_z}}{F_y}(\omega_i) & \frac{D_{1_z}}{F_z}(\omega_i) \\
|
||||
% \frac{D_{2_x}}{F_x}(\omega_i) & \frac{D_{2_x}}{F_y}(\omega_i) & \frac{D_{2_x}}{F_z}(\omega_i) \\
|
||||
% \vdots & \vdots & \vdots \\
|
||||
% \frac{D_{23_z}}{F_x}(\omega_i) & \frac{D_{23_z}}{F_y}(\omega_i) & \frac{D_{23_z}}{F_z}(\omega_i) \\
|
||||
% \end{bmatrix}
|
||||
% \end{equation}
|
||||
% #+end_important
|
||||
|
||||
% We generate such FRF matrix from the measurements using the following script.
|
||||
|
||||
n_meas = 24;
|
||||
n_acc = 23;
|
||||
|
||||
dirs = 'XYZ';
|
||||
|
||||
% Number of Accelerometer * DOF for each acccelerometer / Number of excitation / frequency points
|
||||
FRFs = zeros(3*n_acc, 3, 801);
|
||||
COHs = zeros(3*n_acc, 3, 801);
|
||||
|
||||
% Loop through measurements
|
||||
for i = 1:n_meas
|
||||
% Load the measurement file
|
||||
meas = load(sprintf('mat/meas_frf_coh_%i.mat', i));
|
||||
|
||||
% First: determine what is the exitation (direction and sign)
|
||||
exc_dir = meas.FFT1_AvXSpc_2_1_RMS_RfName(end);
|
||||
exc_sign = meas.FFT1_AvXSpc_2_1_RMS_RfName(end-1);
|
||||
% Determine what is the correct excitation sign
|
||||
exc_factor = str2num([exc_sign, '1']);
|
||||
if exc_dir ~= 'Z'
|
||||
exc_factor = exc_factor*(-1);
|
||||
end
|
||||
|
||||
% Then: loop through the nine measurements and store them at the correct location
|
||||
for j = 2:10
|
||||
% Determine what is the accelerometer and direction
|
||||
[indices_acc_i] = strfind(meas.(sprintf('FFT1_H1_%i_1_RpName', j)), '.');
|
||||
acc_i = str2num(meas.(sprintf('FFT1_H1_%i_1_RpName', j))(indices_acc_i(1)+1:indices_acc_i(2)-1));
|
||||
|
||||
meas_dir = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end);
|
||||
meas_sign = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end-1);
|
||||
% Determine what is the correct measurement sign
|
||||
meas_factor = str2num([meas_sign, '1']);
|
||||
if meas_dir ~= 'Z'
|
||||
meas_factor = meas_factor*(-1);
|
||||
end
|
||||
|
||||
% FRFs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
|
||||
% COHs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
|
||||
|
||||
FRFs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
|
||||
COHs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
|
||||
end
|
||||
end
|
||||
freqs = meas.FFT1_Coh_10_1_RMS_X_Val;
|
||||
|
||||
|
||||
|
||||
% And we save the obtained FRF matrix and Coherence matrix in a =.mat= file.
|
||||
|
||||
save('./mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
|
||||
Reference in New Issue
Block a user