145 lines
5.8 KiB
Matlab
145 lines
5.8 KiB
Matlab
%% Clear Workspace and Close figures
|
|
clear; close all; clc;
|
|
|
|
%% Intialize Laplace variable
|
|
s = zpk('s');
|
|
|
|
% Load the processed data and the model
|
|
|
|
load('./mat/spindle_data.mat', 'spindle');
|
|
load('./mat/spindle_model.mat', 'sys');
|
|
|
|
% Compute the PSD
|
|
|
|
n_av = 4; % Number of average
|
|
|
|
[pxx_1rpm, f_1rpm] = pwelch(spindle.rpm1.xasync, hanning(ceil(length(spindle.rpm1.xasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
[pyy_1rpm, ~] = pwelch(spindle.rpm1.yasync, hanning(ceil(length(spindle.rpm1.yasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
[pzz_1rpm, ~] = pwelch(spindle.rpm1.zasync, hanning(ceil(length(spindle.rpm1.zasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
|
|
[pxx_60rpm, f_60rpm] = pwelch(spindle.rpm60.xasync, hanning(ceil(length(spindle.rpm60.xasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
[pyy_60rpm, ~] = pwelch(spindle.rpm60.yasync, hanning(ceil(length(spindle.rpm60.yasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
[pzz_60rpm, ~] = pwelch(spindle.rpm60.zasync, hanning(ceil(length(spindle.rpm60.zasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
|
|
% Plot the computed PSD
|
|
% The Amplitude Spectral Densities of the displacement of the spindle for the $x$, $y$ and $z$ directions are shown figure [[fig:spindle_psd_xyz_60rpm]]. They correspond to the Asynchronous part shown figure [[fig:spindle_60rpm_sync_async]].
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_1rpm, (pxx_1rpm).^.5, 'DisplayName', '$P_{xx}$ - 1rpm');
|
|
plot(f_1rpm, (pyy_1rpm).^.5, 'DisplayName', '$P_{yy}$ - 1rpm');
|
|
plot(f_1rpm, (pzz_1rpm).^.5, 'DisplayName', '$P_{zz}$ - 1rpm');
|
|
% plot(f_1rpm, spindle.rpm1.adcn*ones(size(f_1rpm)), '--k', 'DisplayName', 'ADC - 1rpm');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
legend('Location', 'northeast');
|
|
xlim([f_1rpm(2), f_1rpm(end)]);
|
|
|
|
|
|
|
|
% #+NAME: fig:spindle_psd_xyz_1rpm
|
|
% #+CAPTION: Power spectral density of the Asynchronous displacement - 1rpm
|
|
% #+RESULTS: fig:spindle_psd_xyz_1rpm
|
|
% [[file:figs/spindle_psd_xyz_1rpm.png]]
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, (pxx_60rpm).^.5, 'DisplayName', '$P_{xx}$ - 60rpm');
|
|
plot(f_60rpm, (pyy_60rpm).^.5, 'DisplayName', '$P_{yy}$ - 60rpm');
|
|
plot(f_60rpm, (pzz_60rpm).^.5, 'DisplayName', '$P_{zz}$ - 60rpm');
|
|
% plot(f_60rpm, spindle.rpm60.adcn*ones(size(f_60rpm)), '--k', 'DisplayName', 'ADC - 60rpm');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
legend('Location', 'northeast');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
|
|
% Compute the response of the model
|
|
|
|
Tfd = abs(squeeze(freqresp(sys('d', 'f'), f_60rpm, 'Hz')));
|
|
|
|
% Plot the PSD of the Force using the model
|
|
|
|
figure;
|
|
plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$P_{xx}$');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
|
|
% Estimated Shape of the PSD of the force
|
|
|
|
s = tf('s');
|
|
|
|
Wd_simple = 1e-8/(1+s/2/pi/0.5)/(1+s/2/pi/100);
|
|
Wf_simple = Wd_simple/tf(sys('d', 'f'));
|
|
TWf_simple = abs(squeeze(freqresp(Wf_simple, f_60rpm, 'Hz')));
|
|
|
|
% Wf = 0.48902*(s+327.9)*(s^2 + 109.6*s + 1.687e04)/((s^2 + 30.59*s + 8541)*(s^2 + 29.11*s + 3.268e04));
|
|
% Wf = 0.15788*(s+418.6)*(s+1697)^2*(s^2 + 124.3*s + 2.529e04)*(s^2 + 681.3*s + 9.018e05)/((s^2 + 23.03*s + 8916)*(s^2 + 33.85*s + 6.559e04)*(s^2 + 71.43*s + 4.283e05)*(s^2 + 40.64*s + 1.789e06));
|
|
|
|
Wf = (s+1697)^2*(s^2 + 114.5*s + 2.278e04)*(s^2 + 205.1*s + 1.627e05)*(s^2 + 285.8*s + 8.624e05)*(s+100)/((s+0.5)*3012*(s^2 + 23.03*s + 8916)*(s^2 + 17.07*s + 4.798e04)*(s^2 + 41.17*s + 4.347e05)*(s^2 + 78.99*s + 1.789e06));
|
|
|
|
TWf = abs(squeeze(freqresp(Wf, f_60rpm, 'Hz')));
|
|
|
|
% PSD in [N]
|
|
% Above $200Hz$, the Amplitude Spectral Density seems dominated by noise coming from the electronics (charge amplifier, ADC, ...). So we don't know what is the frequency content of the force above that frequency. However, we assume that $P_{xx}$ is decreasing with $1/f$ as it seems so be the case below $100Hz$ (figure [[fig:spindle_psd_xyz_60rpm]]).
|
|
|
|
% We then fit the PSD of the displacement with a transfer function (figure [[fig:spindle_psd_d_comp_60rpm]]).
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$\sqrt{P_{xx}}/|T_{d/f}|$');
|
|
plot(f_60rpm, TWf, 'DisplayName', 'Wf');
|
|
plot(f_60rpm, TWf_simple, '-k', 'DisplayName', 'Wfs');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'northeast');
|
|
|
|
% PSD in [m]
|
|
% To obtain the PSD of the force $f$ that induce such displacement, we use the following formula:
|
|
% \[ \sqrt{PSD(d)} = |T_{d/f}| \sqrt{PSD(f)} \]
|
|
|
|
% And so we have:
|
|
% \[ \sqrt{PSD(f)} = |T_{d/f}|^{-1} \sqrt{PSD(d)} \]
|
|
|
|
% The obtain Power Spectral Density of the force is displayed figure [[fig:spindle_psd_f_comp_60rpm]].
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, pxx_60rpm.^.5, 'DisplayName', '$\sqrt{P_{xx}}$');
|
|
plot(f_60rpm, TWf.*Tfd, 'DisplayName', '$|W_f|*|T_{d/f}|$');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'northeast');
|
|
|
|
% Compute the resulting RMS value [m]
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, 1e9*cumtrapz(f_60rpm, (pxx_60rpm)).^.5, '--', 'DisplayName', 'Exp. Data');
|
|
plot(f_60rpm, 1e9*cumtrapz(f_60rpm, ((TWf.*Tfd).^2)).^.5, '--', 'DisplayName', 'Estimated');
|
|
hold off;
|
|
set(gca, 'XScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'southeast');
|
|
|
|
% Compute the resulting RMS value [m]
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f_1rpm, 1e9*cumtrapz(f_1rpm, (pxx_1rpm)), '--', 'DisplayName', 'Exp. Data');
|
|
plot(f_1rpm, 1e9*(f_1rpm(end)-f_1rpm(1))/(length(f_1rpm)-1)*cumsum(pxx_1rpm), '--', 'DisplayName', 'Exp. Data');
|
|
hold off;
|
|
set(gca, 'XScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
|
|
xlim([f_1rpm(2), f_1rpm(end)]);
|
|
legend('Location', 'southeast');
|