diff --git a/figs/frf_meas_sin_excitation.pdf b/figs/frf_meas_sin_excitation.pdf new file mode 100644 index 0000000..88196f1 Binary files /dev/null and b/figs/frf_meas_sin_excitation.pdf differ diff --git a/figs/frf_meas_sin_excitation.png b/figs/frf_meas_sin_excitation.png new file mode 100644 index 0000000..2bb99d1 Binary files /dev/null and b/figs/frf_meas_sin_excitation.png differ diff --git a/matlab/frf_analyze.m b/matlab/frf_analyze.m index c497aac..bf6904e 100644 --- a/matlab/frf_analyze.m +++ b/matlab/frf_analyze.m @@ -1,52 +1,36 @@ -% Analysis -% :PROPERTIES: -% :header-args: :tangle matlab/frf_analyze.m -% :END: +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./src/'); + +% Test with one APA + +%% Load measurement data for APA number 1 +load(sprintf('mat/frf_data_%i.mat', 1), 't', 'Va', 'Vs', 'de', 'da'); -%% Load all the measurements -% meas_data = {}; -% for i = 1:7 -% meas_data(i) = {load(sprintf('mat/frf_data_%i.mat', i), 't', 'Va', 'Vs', 'da', 'de')}; -% end -%% -load(sprintf('mat/frf_data_%i_sweep.mat', 1), 't', 'Va', 'Vs', 'da', 'de') +% Time domain data: -%% figure; plot(t, de); -%% -figure; -plot(t, Va); -%% + +% Compute transfer functions: + Ts = (t(end) - t(1))/(length(t)-1); Fs = 1/Ts; -win = hanning(ceil(5*Fs)); % Hannning Windows +win = hanning(ceil(0.5*Fs)); % Hannning Windows -%% [G_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts); [G_d, ~] = tfestimate(Va, da, win, [], [], 1/Ts); [G_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts); -[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts); -[coh_d, ~] = mscohere(Va, da, win, [], [], 1/Ts); -[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts); - -%% -figure; -hold on; -plot(f, coh_dvf); -plot(f, coh_d); -plot(f, coh_iff); -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlim([1, 5e3]); ylim([0, 1]); - -%% figure; tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); @@ -58,6 +42,7 @@ hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Amplitude $V_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]); hold off; +ylim([10, 30]); ax2 = nexttile; hold on; @@ -72,7 +57,6 @@ yticks(-360:90:360); linkaxes([ax1,ax2],'x'); xlim([5, 5e3]); -%% figure; tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); @@ -81,6 +65,7 @@ plot(f, abs(G_iff)); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Amplitude $V_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]); hold off; +ylim([10, 30]); ax2 = nexttile; plot(f, 180/pi*angle(G_iff)); @@ -90,4 +75,12 @@ hold off; yticks(-360:90:360); linkaxes([ax1,ax2],'x'); -xlim([5, 5e3]); \ No newline at end of file +xlim([5, 5e3]); + +% Comparison of all APA + +%% Load all the measurements +meas_data = {}; +for i = 1:7 + meas_data(i) = {load(sprintf('mat/frf_data_%i.mat', i), 't', 'Va', 'Vs', 'de', 'da')}; +end diff --git a/matlab/frf_save.m b/matlab/frf_save.m index 1bf41ec..33498a1 100644 --- a/matlab/frf_save.m +++ b/matlab/frf_save.m @@ -1,4 +1,4 @@ -% Save Data +% =frf_save.m= - Save Data % :PROPERTIES: % :header-args: :tangle matlab/frf_save.m % :END: diff --git a/matlab/frf_setup.m b/matlab/frf_setup.m index c5072a1..f0e8046 100644 --- a/matlab/frf_setup.m +++ b/matlab/frf_setup.m @@ -1,25 +1,32 @@ -s = tf('s'); -addpath('src') +%% Clear Workspace and Close figures +clear; close all; clc; -%% -Fs = 10e3; % Sampling Frequency [Hz] -Ts = 1/Fs; % Sampling Time [s] +%% Intialize Laplace variable +s = zpk('s'); +addpath('./src/'); + +%% Simulation configuration +Fs = 10e3; % Sampling Frequency [Hz] +Ts = 1/Fs; % Sampling Time [s] + +%% Data record configuration Trec_start = 5; % Start time for Recording [s] Trec_dur = 100; % Recording Duration [s] -Tsim = 2*Trec_start + Trec_dur; % Simulation Time [s] +Tsim = 2*Trec_start + Trec_dur; % Simulation Time [s] %% Sweep Sine gc = 0.1; xi = 0.5; wn = 2*pi*94.3; +% Notch filter at the resonance of the APA G_sweep = 0.2*(s^2 + 2*gc*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2); V_sweep = generateSweepExc('Ts', Ts, ... 'f_start', 10, ... - 'f_end', 2e3, ... + 'f_end', 1e3, ... 'V_mean', 3.25, ... 't_start', Trec_start, ... 'exc_duration', Trec_dur, ... @@ -32,7 +39,16 @@ V_noise = generateShapedNoise('Ts', 1/Fs, ... 't_start', Trec_start, ... 'exc_duration', Trec_dur, ... 'smooth_ends', true, ... - 'V_exc', 0.00/(1 + s/2/pi/50)); + 'V_exc', 0.05/(1 + s/2/pi/10)); + +%% Sinus excitation with increasing amplitude +V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ... + 'V_mean', 3.25, ... + 'sin_ampls', [0.1, 0.2, 0.4, 1, 2, 4], ... + 'sin_period', 1, ... + 'sin_num', 5, ... + 't_start', 10, ... + 'smooth_ends', true); %% Select the excitation signal V_exc = timeseries(V_noise(2,:), V_noise(1,:)); @@ -41,16 +57,16 @@ figure; tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None'); ax1 = nexttile; -plot(V_exc.Time, squeeze(V_exc.Data)); +plot(V_exc(1,:), V_exc(2,:)); xlabel('Time [s]'); ylabel('Amplitude [V]'); ax2 = nexttile; -win = hanning(floor(length(V_exc.Data)/8)); -[pxx, f] = pwelch(squeeze(V_exc.Data), win, 0, [], Fs); +win = hanning(floor(length(V_exc)/8)); +[pxx, f] = pwelch(V_exc(2,:), win, 0, [], Fs); plot(f, pxx) xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]'); set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); xlim([1, Fs/2]); ylim([1e-10, 1e0]); -%% Save +%% Save data that will be loaded in the Simulink file save('./frf_data.mat', 'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', 'V_exc'); diff --git a/matlab/src/generateSinIncreasingAmpl.m b/matlab/src/generateSinIncreasingAmpl.m new file mode 100644 index 0000000..4ce9170 --- /dev/null +++ b/matlab/src/generateSinIncreasingAmpl.m @@ -0,0 +1,54 @@ +function [U_exc] = generateSinIncreasingAmpl(args) +% generateSinIncreasingAmpl - Generate Sinus with increasing amplitude +% +% Syntax: [U_exc] = generateSinIncreasingAmpl(args) +% +% Inputs: +% - args - Optinal arguments: +% - Ts - Sampling Time - [s] +% - V_mean - Mean value of the excitation voltage - [V] +% - sin_ampls - Excitation Amplitudes - [V] +% - sin_freq - Excitation Frequency - [Hz] +% - sin_num - Number of period for each amplitude - [-] +% - t_start - Time at which the excitation begins - [s] +% - smooth_ends - 'true' or 'false': smooth transition between 0 and V_mean - [-] + +arguments + args.Ts (1,1) double {mustBeNumeric, mustBePositive} = 1e-4 + args.V_mean (1,1) double {mustBeNumeric} = 0 + args.sin_ampls double {mustBeNumeric, mustBePositive} = [0.1, 0.2, 0.3] + args.sin_period (1,1) double {mustBeNumeric, mustBePositive} = 1 + args.sin_num (1,1) double {mustBeNumeric, mustBePositive, mustBeInteger} = 3 + args.t_start (1,1) double {mustBeNumeric, mustBePositive} = 5 + args.smooth_ends logical {mustBeNumericOrLogical} = true +end + +t_noise = 0:args.Ts:args.sin_period*args.sin_num; +sin_exc = []; + +for sin_ampl = args.sin_ampls + sin_exc = [sin_exc, args.V_mean + sin_ampl*sin(2*pi/args.sin_period*t_noise)]; +end + +t_smooth_start = args.Ts:args.Ts:args.t_start; + +V_smooth_start = zeros(size(t_smooth_start)); +V_smooth_end = zeros(size(t_smooth_start)); + +if args.smooth_ends + Vd_max = args.V_mean/(0.7*args.t_start); + + V_d = zeros(size(t_smooth_start)); + V_d(t_smooth_start < 0.2*args.t_start) = t_smooth_start(t_smooth_start < 0.2*args.t_start)*Vd_max/(0.2*args.t_start); + V_d(t_smooth_start > 0.2*args.t_start & t_smooth_start < 0.7*args.t_start) = Vd_max; + V_d(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) = Vd_max - (t_smooth_start(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) - 0.7*args.t_start)*Vd_max/(0.2*args.t_start); + + V_smooth_start = cumtrapz(V_d)*args.Ts; + + V_smooth_end = args.V_mean - V_smooth_start; +end + +V_exc = [V_smooth_start, sin_exc, V_smooth_end]; +t_exc = args.Ts*[0:1:length(V_exc)-1]; + +U_exc = [t_exc; V_exc]; diff --git a/test-bench-apa300ml.org b/test-bench-apa300ml.org index ac81355..c227ce6 100644 --- a/test-bench-apa300ml.org +++ b/test-bench-apa300ml.org @@ -25,6 +25,7 @@ #+PROPERTY: header-args:matlab+ :results none #+PROPERTY: header-args:matlab+ :eval no-export #+PROPERTY: header-args:matlab+ :noweb yes +#+PROPERTY: header-args:matlab+ :tangle no #+PROPERTY: header-args:matlab+ :mkdirp yes #+PROPERTY: header-args:matlab+ :output-dir figs @@ -717,18 +718,29 @@ addpath('./matlab/'); addpath('./src/'); #+end_src -** Measurement Setup +** =frf_setup.m= - Measurement Setup :PROPERTIES: :header-args:matlab: :tangle matlab/frf_setup.m :header-args:matlab+: :comments no :END: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) +<> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes +<> +#+end_src + +#+begin_src matlab :eval no :exports none + +addpath('./src/'); +#+end_src First is defined the sampling frequency: #+begin_src matlab %% Simulation configuration Fs = 10e3; % Sampling Frequency [Hz] Ts = 1/Fs; % Sampling Time [s] -Tsim = 110; % Simulation Time [s] #+end_src #+begin_src matlab @@ -737,9 +749,20 @@ Trec_start = 5; % Start time for Recording [s] Trec_dur = 100; % Recording Duration [s] #+end_src +#+begin_src matlab +Tsim = 2*Trec_start + Trec_dur; % Simulation Time [s] +#+end_src + The maximum excitation voltage at resonance is 9Vrms, therefore corresponding to 0.6V of output DAC voltage. #+begin_src matlab %% Sweep Sine +gc = 0.1; +xi = 0.5; +wn = 2*pi*94.3; + +% Notch filter at the resonance of the APA +G_sweep = 0.2*(s^2 + 2*gc*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2); + V_sweep = generateSweepExc('Ts', Ts, ... 'f_start', 10, ... 'f_end', 1e3, ... @@ -747,7 +770,7 @@ V_sweep = generateSweepExc('Ts', Ts, ... 't_start', Trec_start, ... 'exc_duration', Trec_dur, ... 'sweep_type', 'log', ... - 'V_exc', 0.5/(1 + s/2/pi/100)); + 'V_exc', G_sweep*1/(1 + s/2/pi/500)); #+end_src #+begin_src matlab :exports none :tangle no @@ -815,10 +838,36 @@ exportFig('figs/frf_meas_noise_excitation.pdf', 'width', 'full', 'height', 'norm #+RESULTS: [[file:figs/frf_meas_noise_excitation.png]] +#+begin_src matlab +%% Sinus excitation with increasing amplitude +V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ... + 'V_mean', 3.25, ... + 'sin_ampls', [0.1, 0.2, 0.4, 1, 2, 4], ... + 'sin_period', 1, ... + 'sin_num', 5, ... + 't_start', 10, ... + 'smooth_ends', true); +#+end_src + +#+begin_src matlab :exports none :tangle no +figure; +plot(V_sin(1,:), V_sin(2,:)); +xlabel('Time [s]'); ylabel('Amplitude [V]'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace +exportFig('figs/frf_meas_sin_excitation.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:frf_meas_sin_excitation +#+caption: Example of Shaped noise excitation signal +#+RESULTS: +[[file:figs/frf_meas_sin_excitation.png]] + Then, we select the wanted excitation signal. #+begin_src matlab %% Select the excitation signal -V_exc = V_noise; +V_exc = timeseries(V_noise(2,:), V_noise(1,:)); #+end_src #+begin_src matlab :exports none :eval no @@ -838,7 +887,12 @@ set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); xlim([1, Fs/2]); ylim([1e-10, 1e0]); #+end_src -** Save Data +#+begin_src matlab +%% Save data that will be loaded in the Simulink file +save('./frf_data.mat', 'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', 'V_exc'); +#+end_src + +** =frf_save.m= - Save Data :PROPERTIES: :header-args: :tangle matlab/frf_save.m :END: @@ -856,10 +910,10 @@ And we load the data on the Workspace: #+begin_src matlab data = SimulinkRealTime.utils.getFileScopeData('data/data.dat').data; -Va = data(:, 1); % Excitation Voltage (input of PD200) [V] -Vs = data(:, 2); % Measured voltage (force sensor) [V] -de = data(:, 3); % Measurment displacement (encoder) [m] -da = data(:, 4); % Measurement displacement (attocube) [m] +da = data(:, 1); % Excitation Voltage (input of PD200) [V] +de = data(:, 2); % Measured voltage (force sensor) [V] +Vs = data(:, 3); % Measurment displacement (encoder) [m] +Va = data(:, 4); % Measurement displacement (attocube) [m] t = data(:, end); % Time [s] #+end_src @@ -867,10 +921,347 @@ And we save this to a =mat= file: #+begin_src matlab apa_number = 1; -save(sprintf('mat/frf_data_%i.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da'); +save(sprintf('mat/frf_data_%i_huddle.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da'); #+end_src -** Analysis +** Measurements on APA 1 +*** Introduction :ignore: +Measurements are first performed on the APA number 1. + +*** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) +<> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes +<> +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab/mat/'); +addpath('./matlab/src/'); +addpath('./matlab/'); +#+end_src + +#+begin_src matlab :eval no +addpath('./mat/'); +addpath('./src/'); +#+end_src + +*** Huddle Test +#+begin_src matlab +load(sprintf('frf_data_%i_huddle.mat', 1), 't', 'Va', 'Vs', 'de', 'da'); +#+end_src + +#+begin_src matlab :exports none +figure; +tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(t, da, 'DisplayName', '$d_a$ - Attocube') +plot(t, de, 'DisplayName', '$d_e$ - Encoder') +hold off; +xlabel('Time [s]'); +ylabel('Displacement [m]'); +legend('location', 'northeast') + +ax2 = nexttile; +hold on; +plot(t, Vs, 'DisplayName', '$V_s$ - Force Sensor') +hold off; +xlabel('Time [s]'); +ylabel('Voltage [V]'); +legend('location', 'northeast') +#+end_src + +#+begin_src matlab +%% Parameter for Spectral analysis +Ts = (t(end) - t(1))/(length(t)-1); +Fs = 1/Ts; + +win = hanning(ceil(10*Fs)); % Hannning Windows +#+end_src + +#+begin_src matlab +[phh_da, f] = pwelch(da - mean(da), win, [], [], 1/Ts); +[phh_de, ~] = pwelch(de - mean(de), win, [], [], 1/Ts); +[phh_Vs, ~] = pwelch(Vs - mean(Vs), win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +figure; +tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(f, phh_da, 'DisplayName', '$d_a$ - Attocube') +plot(f, phh_de, 'DisplayName', '$d_e$ - Encoder') +hold off; +xlabel('Time [s]'); ylabel('ASD [$m/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') +xlim([5e-1, 5e3]); + +ax2 = nexttile; +hold on; +plot(f, phh_Vs, 'DisplayName', '$V_s$ - Force Sensor') +hold off; +xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') +xlim([5e-1, 5e3]); +#+end_src + +*** First identification with Noise +#+begin_src matlab +load(sprintf('mat/frf_data_%i_noise.mat', 1), 't', 'Va', 'Vs', 'da', 'de') +#+end_src + +#+begin_src matlab +[pxx_da, f] = pwelch(da, win, [], [], 1/Ts); +[pxx_de, ~] = pwelch(de, win, [], [], 1/Ts); +[pxx_Vs, ~] = pwelch(Vs, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +figure; +tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(f, sqrt(phh_da), 'DisplayName', 'Huddle') +plot(f, sqrt(pxx_da), 'DisplayName', 'Noise') +hold off; +xlabel('Time [s]'); ylabel('ASD Attocube [$m/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') + +ax2 = nexttile; +hold on; +plot(f, sqrt(phh_Vs), 'DisplayName', 'Huddle') +plot(f, sqrt(pxx_Vs), 'DisplayName', 'Noise') +hold off; +xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') + +linkaxes([ax1,ax2],'x'); +xlim([5e-1, 5e3]); +#+end_src + +#+begin_src matlab +[G_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts); +[G_d, ~] = tfestimate(Va, da, win, [], [], 1/Ts); +[G_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts); + +[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts); +[coh_d, ~] = mscohere(Va, da, win, [], [], 1/Ts); +[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +hold on; +plot(f, coh_dvf); +plot(f, coh_d); +plot(f, coh_iff); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlim([1, 5e3]); ylim([0, 1]); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(f, abs(G_d), 'DisplayName', 'Attocube'); +plot(f, abs(G_dvf), 'DisplayName', 'Encoder'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-8, 1e-3]); +legend('location', 'northeast'); + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_d)); +plot(f, 180/pi*angle(G_dvf)); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); + +linkaxes([ax1,ax2],'x'); +xlim([5e-1, 5e3]); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +plot(f, abs(G_iff)); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); +hold off; + +ax2 = nexttile; +plot(f, 180/pi*angle(G_iff)); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); + +linkaxes([ax1,ax2],'x'); +xlim([1, 2e3]); +#+end_src + +*** Second identification with Sweep +#+begin_src matlab +load(sprintf('mat/frf_data_%i_sweep.mat', 1), 't', 'Va', 'Vs', 'da', 'de') +#+end_src + +#+begin_src matlab +[pxx_da, f] = pwelch(da, win, [], [], 1/Ts); +[pxx_de, ~] = pwelch(de, win, [], [], 1/Ts); +[pxx_Vs, ~] = pwelch(Vs, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +figure; +tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(f, sqrt(phh_da), 'DisplayName', 'Huddle') +plot(f, sqrt(pxx_da), 'DisplayName', 'Noise') +hold off; +xlabel('Time [s]'); ylabel('ASD Attocube [$m/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') + +ax2 = nexttile; +hold on; +plot(f, sqrt(phh_Vs), 'DisplayName', 'Huddle') +plot(f, sqrt(pxx_Vs), 'DisplayName', 'Noise') +hold off; +xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +legend('location', 'northeast') + +linkaxes([ax1,ax2],'x'); +xlim([1e1, 2e3]); +#+end_src + +#+begin_src matlab +[G_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts); +[G_d, ~] = tfestimate(Va, da, win, [], [], 1/Ts); +[G_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts); + +[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts); +[coh_d, ~] = mscohere(Va, da, win, [], [], 1/Ts); +[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +hold on; +plot(f, coh_dvf); +plot(f, coh_d); +plot(f, coh_iff); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlim([10, 2e3]); ylim([0, 1]); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +hold on; +plot(f, abs(G_d), 'DisplayName', 'Attocube'); +plot(f, abs(G_dvf), 'DisplayName', 'Encoder'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-8, 1e-3]); +legend('location', 'northeast'); + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_d)); +plot(f, 180/pi*angle(G_dvf)); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); + +linkaxes([ax1,ax2],'x'); +xlim([10, 2e3]); +#+end_src + +#+begin_src matlab :exports none +%% +figure; +tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile; +plot(f, abs(G_iff)); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); +hold off; + +ax2 = nexttile; +plot(f, 180/pi*angle(G_iff)); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); + +linkaxes([ax1,ax2],'x'); +xlim([10, 2e3]); +#+end_src + +*** Extract Parameters + +Quasi static gain between $d$ and $V_a$: +#+begin_src matlab +g_d_Va = mean(abs(G_dvf(f > 10 & f < 15))); +#+end_src + +#+begin_src matlab :results value replace :exports results +sprintf('g_d_Va = %.1e [m/V]', g_d_Va) +#+end_src + +#+RESULTS: +: g_d_Va = 1.7e-05 [m/V] + +Quasi static gain between $V_s$ and $V_a$: +#+begin_src matlab +g_Vs_Va = mean(abs(G_iff(f > 10 & f < 15))); +#+end_src + +#+begin_src matlab :results value replace :exports results +sprintf('g_Vs_Va = %.1e [V/V]', g_Vs_Va) +#+end_src + +#+RESULTS: +: g_Vs_Va = 5.7e-01 [V/V] + +** Comparison of all APA :PROPERTIES: :header-args: :tangle matlab/frf_analyze.m :END: @@ -3038,5 +3429,106 @@ t_exc = args.Ts*[0:1:length(V_exc)-1]; U_exc = [t_exc; V_exc]; #+end_src +** =generateSinIncreasingAmpl=: Generate Sinus with increasing amplitude +:PROPERTIES: +:header-args:matlab+: :tangle ./matlab/src/generateSinIncreasingAmpl.m +:header-args:matlab+: :comments none :mkdirp yes :eval no +:END: +<> + +*** Function description +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +function [U_exc] = generateSinIncreasingAmpl(args) +% generateSinIncreasingAmpl - Generate Sinus with increasing amplitude +% +% Syntax: [U_exc] = generateSinIncreasingAmpl(args) +% +% Inputs: +% - args - Optinal arguments: +% - Ts - Sampling Time - [s] +% - V_mean - Mean value of the excitation voltage - [V] +% - sin_ampls - Excitation Amplitudes - [V] +% - sin_freq - Excitation Frequency - [Hz] +% - sin_num - Number of period for each amplitude - [-] +% - t_start - Time at which the excitation begins - [s] +% - smooth_ends - 'true' or 'false': smooth transition between 0 and V_mean - [-] +#+end_src + +*** Optional Parameters +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +arguments + args.Ts (1,1) double {mustBeNumeric, mustBePositive} = 1e-4 + args.V_mean (1,1) double {mustBeNumeric} = 0 + args.sin_ampls double {mustBeNumeric, mustBePositive} = [0.1, 0.2, 0.3] + args.sin_period (1,1) double {mustBeNumeric, mustBePositive} = 1 + args.sin_num (1,1) double {mustBeNumeric, mustBePositive, mustBeInteger} = 3 + args.t_start (1,1) double {mustBeNumeric, mustBePositive} = 5 + args.smooth_ends logical {mustBeNumericOrLogical} = true +end +#+end_src + +*** Sinus excitation +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +t_noise = 0:args.Ts:args.sin_period*args.sin_num; +sin_exc = []; +#+end_src + +#+begin_src matlab +for sin_ampl = args.sin_ampls + sin_exc = [sin_exc, args.V_mean + sin_ampl*sin(2*pi/args.sin_period*t_noise)]; +end +#+end_src + +*** Smooth Ends +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +t_smooth_start = args.Ts:args.Ts:args.t_start; + +V_smooth_start = zeros(size(t_smooth_start)); +V_smooth_end = zeros(size(t_smooth_start)); + +if args.smooth_ends + Vd_max = args.V_mean/(0.7*args.t_start); + + V_d = zeros(size(t_smooth_start)); + V_d(t_smooth_start < 0.2*args.t_start) = t_smooth_start(t_smooth_start < 0.2*args.t_start)*Vd_max/(0.2*args.t_start); + V_d(t_smooth_start > 0.2*args.t_start & t_smooth_start < 0.7*args.t_start) = Vd_max; + V_d(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) = Vd_max - (t_smooth_start(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) - 0.7*args.t_start)*Vd_max/(0.2*args.t_start); + + V_smooth_start = cumtrapz(V_d)*args.Ts; + + V_smooth_end = args.V_mean - V_smooth_start; +end +#+end_src + +*** Combine Excitation signals +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +V_exc = [V_smooth_start, sin_exc, V_smooth_end]; +t_exc = args.Ts*[0:1:length(V_exc)-1]; +#+end_src + +#+begin_src matlab +U_exc = [t_exc; V_exc]; +#+end_src + * Bibliography :ignore: #+latex: \printbibliography