Add measurement of spurious resonances

piezostack_apa300_mel_lab_5may2021/piezo_stack_mel_5may2021.m

@@ -0,0 +1,23 @@
+% Vibration measurement on:
+% Piezo stack in MEL lab - 5 may 2021
+
+% Excitation: hammer dytran with metal tip : 230e-6V/N
+% Response: laser vibro dual fibres 1mm/s/V
+
+% recorded signals and saved traces
+% laser fibres on both sides edges of stack 
+--> photo: piezo_rotation.jpg
+--> plot: piezo_stack_rotation_5may2021.png
+% meas3 hammer on edge right: green
+% meas4 hammer on centre top : blue
+% meas5 hammer on edge left : orange
+
+% laser fibres on base and top of stack (vertical axis)
+--> photo: piezo_vertical.jpg
+--> plot: piezo_stack_rotation_vertical_5may2021.png
+% meas7 hammer on vertical centre top : orange
+
+% Stack rotated 90deg - laser fibres on upper interface of stack (top and bottom of interface plate, ~1cm apart)
+--> photo: none
+--> plot: piezo_stack_rotation_interface_5may2021.png
+% meas8 hammer on interface top : orange

test-bench-apa300ml.org

@@ -230,79 +230,6 @@ The excitation frequency is set to be 1kHz.
 There is clearly a problem with APA300ML number 3
 #+end_warning
 
-* Stiffness measurement
-** Matlab Init                                              :noexport:ignore:
-#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
-<<matlab-dir>>
-#+end_src
-
-#+begin_src matlab :exports none :results silent :noweb yes
-<<matlab-init>>
-#+end_src
-
-#+begin_src matlab :tangle no
-addpath('./matlab/mat/');
-addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab :eval no
-addpath('./mat/');
-#+end_src
-
-** APA test
-#+begin_src matlab
-load('meas_stiff_apa_1_x.mat', 't', 'F', 'd');
-#+end_src
-
-#+begin_src matlab
-figure;
-plot(t, F)
-#+end_src
-
-#+begin_src matlab
-%% Automatic Zero of the force
-F = F - mean(F(t > 0.1 & t < 0.3));
-
-%% Start measurement at t = 0.2 s
-d = d(t > 0.2);
-F = F(t > 0.2);
-t = t(t > 0.2); t = t - t(1);
-#+end_src
-
-#+begin_src matlab
-i_l_start = find(F > 0.3, 1, 'first');
-[~, i_l_stop] = max(F);
-#+end_src
-
-#+begin_src matlab
-F_l = F(i_l_start:i_l_stop);
-d_l = d(i_l_start:i_l_stop);
-#+end_src
-
-#+begin_src matlab
-fit_l = polyfit(F_l, d_l, 1);
-
-% %% Reset displacement based on fit
-% d = d - fit_l(2);
-% fit_s(2) = fit_s(2) - fit_l(2);
-% fit_l(2) = 0;
-
-% %% Estimated Stroke
-% F_max = fit_s(2)/(fit_l(1) - fit_s(1));
-% d_max = fit_l(1)*F_max;
-#+end_src
-
-#+begin_src matlab
-h^2/fit_l(1)
-#+end_src
-
-#+begin_src matlab
-figure;
-hold on;
-plot(F,d,'k')
-plot(F_l, d_l)
-plot(F_l, F_l*fit_l(1) + fit_l(2), '--')
-#+end_src
 * Stroke measurement
 ** Introduction                                                      :ignore:
 
@@ -564,6 +491,79 @@ data2orgtable(1e6*apa300ml_stroke', {'APA 1', 'APA 2', 'APA 3', 'APA 4', 'APA 5'
 | APA 6 |              363.9 |
 | APA 7 |              358.4 |
 
+* TODO Stiffness measurement
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('./matlab/mat/');
+addpath('./matlab/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('./mat/');
+#+end_src
+
+** APA test
+#+begin_src matlab
+load('meas_stiff_apa_1_x.mat', 't', 'F', 'd');
+#+end_src
+
+#+begin_src matlab
+figure;
+plot(t, F)
+#+end_src
+
+#+begin_src matlab
+%% Automatic Zero of the force
+F = F - mean(F(t > 0.1 & t < 0.3));
+
+%% Start measurement at t = 0.2 s
+d = d(t > 0.2);
+F = F(t > 0.2);
+t = t(t > 0.2); t = t - t(1);
+#+end_src
+
+#+begin_src matlab
+i_l_start = find(F > 0.3, 1, 'first');
+[~, i_l_stop] = max(F);
+#+end_src
+
+#+begin_src matlab
+F_l = F(i_l_start:i_l_stop);
+d_l = d(i_l_start:i_l_stop);
+#+end_src
+
+#+begin_src matlab
+fit_l = polyfit(F_l, d_l, 1);
+
+% %% Reset displacement based on fit
+% d = d - fit_l(2);
+% fit_s(2) = fit_s(2) - fit_l(2);
+% fit_l(2) = 0;
+
+% %% Estimated Stroke
+% F_max = fit_s(2)/(fit_l(1) - fit_s(1));
+% d_max = fit_l(1)*F_max;
+#+end_src
+
+#+begin_src matlab
+h^2/fit_l(1)
+#+end_src
+
+#+begin_src matlab
+figure;
+hold on;
+plot(F,d,'k')
+plot(F_l, d_l)
+plot(F_l, F_l*fit_l(1) + fit_l(2), '--')
+#+end_src
 * Test-Bench Description
 
 #+begin_note
@@ -580,6 +580,7 @@ Here are the documentation of the equipment used for this test bench:
 [[file:figs/test_bench_apa_alone.png]]
 
 * Measurement Procedure
+<<sec:measurement_procedure>>
 ** Introduction                                                      :ignore:
 
 ** Stroke Measurement
@@ -621,7 +622,7 @@ For each excitation amplitude, the vertical displacement $d$ of the mass is meas
 Then, $d$ is plotted as a function of $V_a$ for all the amplitudes.
 
 #+name: fig:expected_hysteresis
-#+caption: Expected Hysteresis (cite:poel10_explor_activ_hard_mount_vibrat)
+#+caption: Expected Hysteresis cite:poel10_explor_activ_hard_mount_vibrat
 #+attr_latex: :width 0.8\linewidth
 [[file:figs/expected_hysteresis.png]]
 
@@ -685,6 +686,298 @@ This can also be performed with and without the encoder fixed to the APA.
 
 Compare all the obtained parameters for all the test APA.
 
+* FRF measurement
+** Introduction                                                      :ignore:
+
+- [ ] Schematic of the measurement
+
+| Variable |                              | Unit | Hardware                      |
+|----------+------------------------------+------+-------------------------------|
+| =Va=     | Output DAC voltage           | [V]  | DAC - Ch. 1 => PD200 => APA   |
+| =Vs=     | Measured stack voltage (ADC) | [V]  | APA => ADC - Ch. 1            |
+| =de=     | Encoder Measurement          | [m]  | PEPU Ch. 1 - IO318(1) - Ch. 1 |
+| =da=     | Attocube Measurement         | [m]  | PEPU Ch. 2 - IO318(1) - Ch. 2 |
+| =t=      | Time                         | [s]  |                               |
+
+** Matlab Init                                             :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('./matlab/src/');
+addpath('./matlab/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('./src/');
+#+end_src
+
+** Measurement Setup
+:PROPERTIES:
+:header-args:matlab: :tangle matlab/frf_setup.m
+:header-args:matlab+: :comments no
+:END:
+
+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
+%% Data record configuration
+Trec_start = 5;  % Start time for Recording [s]
+Trec_dur   = 100; % Recording Duration [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
+V_sweep = generateSweepExc('Ts',      Ts, ...
+                           'f_start', 10, ...
+                           'f_end',   1e3, ...
+                           'V_mean',  3.25, ...
+                           't_start', Trec_start, ...
+                           'exc_duration', Trec_dur, ...
+                           'sweep_type',   'log', ...
+                           'V_exc', 0.5/(1 + s/2/pi/100));
+#+end_src
+
+#+begin_src matlab :exports none :tangle no
+figure;
+tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(V_sweep(1,:), V_sweep(2,:));
+xlabel('Time [s]'); ylabel('Amplitude [V]');
+
+ax2 = nexttile;
+win = hanning(floor(length(V_sweep(2,:))/80));
+[pxx, f] = pwelch(V_sweep(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]);
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/frf_meas_sweep_excitation.pdf', 'width', 'full', 'height', 'normal');
+#+end_src
+
+#+name: fig:frf_meas_sweep_excitation
+#+caption: Example of Sweep Sin excitation signal
+#+RESULTS:
+[[file:figs/frf_meas_sweep_excitation.png]]
+
+
+A white noise excitation signal can be very useful in order to obtain a first idea of the plant FRF.
+The gain can be gradually increased until satisfactory output is obtained.
+#+begin_src matlab
+%% Shaped Noise
+V_noise = generateShapedNoise('Ts', 1/Fs, ...
+                              'V_mean', 3.25, ...
+                              't_start', Trec_start, ...
+                              'exc_duration', Trec_dur, ...
+                              'smooth_ends', true, ...
+                              'V_exc', 0.05/(1 + s/2/pi/10));
+#+end_src
+
+#+begin_src matlab :exports none :tangle no
+figure;
+tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(V_noise(1,:), V_noise(2,:));
+xlabel('Time [s]'); ylabel('Amplitude [V]');
+
+ax2 = nexttile;
+win = hanning(floor(length(V_noise)/8));
+[pxx, f] = pwelch(V_noise(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]);
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/frf_meas_noise_excitation.pdf', 'width', 'full', 'height', 'normal');
+#+end_src
+
+#+name: fig:frf_meas_noise_excitation
+#+caption: Example of Shaped noise excitation signal
+#+RESULTS:
+[[file:figs/frf_meas_noise_excitation.png]]
+
+Then, we select the wanted excitation signal.
+#+begin_src matlab
+%% Select the excitation signal
+V_exc = V_noise;
+#+end_src
+
+#+begin_src matlab :exports none :eval no
+figure;
+tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(V_exc(1,:), V_exc(2,:));
+xlabel('Time [s]'); ylabel('Amplitude [V]');
+
+ax2 = nexttile;
+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]);
+#+end_src
+
+** Save Data
+:PROPERTIES:
+:header-args: :tangle matlab/frf_save.m
+:END:
+
+First, we get data from the Speedgoat:
+#+begin_src matlab
+tg = slrt;
+
+f = SimulinkRealTime.openFTP(tg);
+mget(f, 'data/data.dat');
+close(f);
+#+end_src
+
+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]
+t  = data(:, end); % Time [s]
+#+end_src
+
+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');
+#+end_src
+
+** Analysis
+:PROPERTIES:
+:header-args: :tangle matlab/frf_analyze.m
+:END:
+
+*** Matlab Init                                            :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('./matlab/src/');
+addpath('./matlab/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('./src/');
+#+end_src
+
+*** Test with one APA
+#+begin_src matlab
+%% Load measurement data for APA number 1
+load(sprintf('mat/frf_data_%i.mat', 1), 't', 'Va', 'Vs', 'de', 'da');
+#+end_src
+
+Time domain data:
+#+begin_src matlab
+figure;
+plot(t, de);
+#+end_src
+
+Compute transfer functions:
+#+begin_src matlab
+Ts = (t(end) - t(1))/(length(t)-1);
+Fs = 1/Ts;
+
+win = hanning(ceil(0.5*Fs)); % Hannning Windows
+#+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);
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+hold on;
+plot(f, abs(G_dvf));
+plot(f, abs(G_d));
+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;
+plot(f, 180/pi*angle(G_dvf));
+plot(f, 180/pi*angle(G_d));
+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([5, 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_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+ylim([10, 30]);
+
+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([5, 5e3]);
+#+end_src
+
+*** Comparison of all APA
+#+begin_src matlab :exports none
+%% 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
+#+end_src
+
 * Measurement Results
 
 * TODO Compare with the FEM/Simscape Model                          :noexport:
@@ -2307,7 +2600,7 @@ The APA is excited with an instrumented hammer and the transfer function from th
 #+name: fig:measurement_setup_torsion
 #+caption: Measurement setup with a Laser Doppler Vibrometer and one instrumental hammer
 #+attr_latex: :width 0.7\linewidth
-[[file:figs/measurement_setup_torsion.png]]
+[[file:figs/measurement_setup_torsion.jpg]]
 
 ** Bending - X
 
@@ -2318,20 +2611,21 @@ The impact point is on the back-side of the APA aligned with the top measurement
 #+name: fig:measurement_setup_X_bending
 #+caption: X-Bending measurement setup
 #+attr_latex: :width 0.7\linewidth
-[[file:figs/measurement_setup_X_bending.png]]
+[[file:figs/measurement_setup_X_bending.jpg]]
 
 The data is loaded.
 #+begin_src matlab
-bending_X = load('apa300ml_bending_X_top.mat')
+bending_X = load('apa300ml_bending_X_top.mat');
 #+end_src
 
+The config for =tfestimate= is performed:
 #+begin_src matlab
 Ts = bending_X.Track1_X_Resolution; % Sampling frequency [Hz]
+win = hann(ceil(1/Ts));
 #+end_src
 
 The transfer function from the input force to the output "rotation" (difference between the two measured distances).
 #+begin_src matlab
-win = hann(ceil(1/Ts));
 [G_bending_X, f]  = tfestimate(bending_X.Track1, bending_X.Track2, win, [], [], 1/Ts);
 #+end_src
 
@@ -2369,11 +2663,11 @@ The impact point of the instrumented hammer is located on the back surface of th
 #+name: fig:measurement_setup_Y_bending
 #+caption: Y-Bending measurement setup
 #+attr_latex: :width 0.7\linewidth
-[[file:figs/measurement_setup_Y_bending.png]]
+[[file:figs/measurement_setup_Y_bending.jpg]]
 
 The data is loaded, and the transfer function from the force to the measured rotation is computed.
 #+begin_src matlab
-bending_Y = load('apa300ml_bending_Y_top.mat')
+bending_Y = load('apa300ml_bending_Y_top.mat');
 [G_bending_Y, ~]  = tfestimate(bending_Y.Track1, bending_Y.Track2, win, [], [], 1/Ts);
 #+end_src
 
@@ -2410,12 +2704,12 @@ The excitation is shown on the other side of the APA, on the side to excite the
 #+name: fig:measurement_setup_torsion_bis
 #+caption: Z-Torsion measurement setup
 #+attr_latex: :width 0.7\linewidth
-[[file:figs/measurement_setup_torsion_bis.png]]
+[[file:figs/measurement_setup_torsion_bis.jpg]]
 
 The data is loaded, and the transfer function computed.
 #+begin_src matlab
-torsion = load('apa300ml_torsion_left.mat')
-[G_torsion, ~]  = tfestimate(torsion_left.Track1, torsion_left.Track2, win, [], [], 1/Ts);
+torsion = load('apa300ml_torsion_left.mat');
+[G_torsion, ~]  = tfestimate(torsion.Track1, torsion.Track2, win, [], [], 1/Ts);
 #+end_src
 
 The results are shown in Figure [[fig:apa300ml_meas_freq_torsion_z]].
@@ -2431,7 +2725,8 @@ set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
 xlabel('Frequency [Hz]'); ylabel('Amplitude');
 xlim([50, 2e3]); ylim([1e-5, 2e-2])
 text(415, 4.3e-3,{'415Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
-text(267, 8e-4,{'267Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(267, 8e-4,{'267Hz'}, 'VerticalAlignment', 'bottom','HorizontalAlignment','center')
+text(800, 6e-4,{'800Hz'}, 'VerticalAlignment', 'bottom','HorizontalAlignment','center')
 #+end_src
 
 #+begin_src matlab :tangle no :exports results :results file replace
@@ -2443,6 +2738,39 @@ exportFig('figs/apa300ml_meas_freq_torsion_z.pdf', 'width', 'wide', 'height', 'n
 #+RESULTS:
 [[file:figs/apa300ml_meas_freq_torsion_z.png]]
 
+In order to verify that, the APA is excited on the top part such that the torsion mode should not be excited.
+#+begin_src matlab
+torsion = load('apa300ml_torsion_top.mat');
+[G_torsion_top, ~]  = tfestimate(torsion.Track1, torsion.Track2, win, [], [], 1/Ts);
+#+end_src
+
+The two FRF are compared in Figure [[fig:apa300ml_meas_freq_torsion_z_comp]].
+It is clear that the first two modes does not correspond to the torsional mode.
+Maybe the resonance at 800Hz, or even higher resonances. It is difficult to conclude here.
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, abs(G_torsion), 'k-', 'DisplayName', 'Left excitation');
+plot(f, abs(G_torsion_top), '-', 'DisplayName', 'Top excitation');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([50, 2e3]); ylim([1e-5, 2e-2])
+text(415, 4.3e-3,{'415Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(267, 8e-4,{'267Hz'}, 'VerticalAlignment', 'bottom','HorizontalAlignment','center')
+text(800, 2e-3,{'800Hz'}, 'VerticalAlignment', 'bottom','HorizontalAlignment','center')
+legend('location', 'northwest');
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/apa300ml_meas_freq_torsion_z_comp.pdf', 'width', 'wide', 'height', 'normal');
+#+end_src
+
+#+name: fig:apa300ml_meas_freq_torsion_z_comp
+#+caption: Obtained FRF for the Z-torsion
+#+RESULTS:
+[[file:figs/apa300ml_meas_freq_torsion_z_comp.png]]
+
 ** Compare
 The three measurements are shown in Figure [[fig:apa300ml_meas_freq_compare]].
 #+begin_src matlab :exports none
@@ -2477,7 +2805,7 @@ It is however quite interesting that there is a factor $\approx \sqrt{2}$ betwee
 
 #+name: tab:apa300ml_measured_modes_freq
 #+caption: Measured frequency of the modes
-#+attr_latex: :environment tabularx :width \linewidth :align lX
+#+attr_latex: :environment tabularx :width 0.6\linewidth :align ccc
 #+attr_latex: :center t :booktabs t :float t
 | Mode      | Strut Mode | Measured Frequency |
 |-----------+------------+--------------------|
@@ -2485,5 +2813,230 @@ It is however quite interesting that there is a factor $\approx \sqrt{2}$ betwee
 | Y-Bending | 285Hz      | 410Hz              |
 | Z-Torsion | 400Hz      | ?                  |
 
+* Function
+** =generateSweepExc=: Generate sweep sinus excitation
+:PROPERTIES:
+:header-args:matlab+: :tangle ./matlab/src/generateSweepExc.m
+:header-args:matlab+: :comments none :mkdirp yes :eval no
+:END:
+<<sec:generateSweepExc>>
+
+*** Function description
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+function [U_exc] = generateSweepExc(args)
+% generateSweepExc - Generate a Sweep Sine excitation signal
+%
+% Syntax: [U_exc] = generateSweepExc(args)
+%
+% Inputs:
+%    - args - Optinal arguments:
+%        - Ts              - Sampling Time                                              - [s]
+%        - f_start         - Start frequency of the sweep                               - [Hz]
+%        - f_end           - End frequency of the sweep                                 - [Hz]
+%        - V_mean          - Mean value of the excitation voltage                       - [V]
+%        - V_exc           - Excitation Amplitude for the Sweep, could be numeric or TF - [V]
+%        - t_start         - Time at which the sweep begins                             - [s]
+%        - exc_duration    - Duration of the sweep                                      - [s]
+%        - sweep_type      - 'logarithmic' or 'linear'                                  - [-]
+%        - 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.f_start         (1,1) double  {mustBeNumeric, mustBePositive} = 1
+    args.f_end           (1,1) double  {mustBeNumeric, mustBePositive} = 1e3
+    args.V_mean          (1,1) double  {mustBeNumeric} = 0
+    args.V_exc                                         = 1
+    args.t_start         (1,1) double  {mustBeNumeric, mustBeNonnegative} = 5
+    args.exc_duration    (1,1) double  {mustBeNumeric, mustBePositive} = 10
+    args.sweep_type            char    {mustBeMember(args.sweep_type,{'log', 'lin'})} = 'lin'
+    args.smooth_ends           logical {mustBeNumericOrLogical} = true
+end
+#+end_src
+
+*** Sweep Sine part
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+t_sweep = 0:args.Ts:args.exc_duration;
+
+if strcmp(args.sweep_type, 'log')
+    V_exc = sin(2*pi*args.f_start * args.exc_duration/log(args.f_end/args.f_start) * (exp(log(args.f_end/args.f_start)*t_sweep/args.exc_duration) - 1));
+elseif strcmp(args.sweep_type, 'lin')
+    V_exc = sin(2*pi*(args.f_start + (args.f_end - args.f_start)/2/args.exc_duration*t_sweep).*t_sweep);
+else
+    error('sweep_type should either be equal to "log" or to "lin"');
+end
+#+end_src
+
+#+begin_src matlab
+if isnumeric(args.V_exc)
+    V_sweep = args.V_mean + args.V_exc*V_exc;
+elseif isct(args.V_exc)
+    if strcmp(args.sweep_type, 'log')
+        V_sweep = args.V_mean + abs(squeeze(freqresp(args.V_exc, args.f_start*(args.f_end/args.f_start).^(t_sweep/args.exc_duration), 'Hz')))'.*V_exc;
+    elseif strcmp(args.sweep_type, 'lin')
+        V_sweep = args.V_mean + abs(squeeze(freqresp(args.V_exc, args.f_start+(args.f_end-args.f_start)/args.exc_duration*t_sweep, 'Hz')))'.*V_exc;
+    end
+end
+#+end_src
+
+*** Smooth Ends
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+if args.t_start > 0
+    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
+else
+    V_smooth_start = [];
+    V_smooth_end = [];
+end
+#+end_src
+
+*** Combine Excitation signals
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+V_exc = [V_smooth_start, V_sweep, 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
+
+** =generateShapedNoise=: Generate Shaped Noise excitation
+:PROPERTIES:
+:header-args:matlab+: :tangle ./matlab/src/generateShapedNoise.m
+:header-args:matlab+: :comments none :mkdirp yes :eval no
+:END:
+<<sec:generateShapedNoise>>
+
+*** Function description
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+function [U_exc] = generateShapedNoise(args)
+% generateShapedNoise - Generate a Shaped Noise excitation signal
+%
+% Syntax: [U_exc] = generateShapedNoise(args)
+%
+% Inputs:
+%    - args - Optinal arguments:
+%        - Ts              - Sampling Time                                              - [s]
+%        - V_mean          - Mean value of the excitation voltage                       - [V]
+%        - V_exc           - Excitation Amplitude, could be numeric or TF               - [V rms]
+%        - t_start         - Time at which the noise begins                             - [s]
+%        - exc_duration    - Duration of the noise                                      - [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.V_exc                                         = 1
+    args.t_start         (1,1) double  {mustBeNumeric, mustBePositive} = 5
+    args.exc_duration    (1,1) double  {mustBeNumeric, mustBePositive} = 10
+    args.smooth_ends           logical {mustBeNumericOrLogical} = true
+end
+#+end_src
+
+*** Shaped Noise
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+t_noise = 0:args.Ts:args.exc_duration;
+
+#+end_src
+
+#+begin_src matlab
+if isnumeric(args.V_exc)
+    V_noise = args.V_mean + args.V_exc*sqrt(1/args.Ts/2)*randn(length(t_noise), 1)';
+elseif isct(args.V_exc)
+    V_noise = args.V_mean + lsim(args.V_exc, sqrt(1/args.Ts/2)*randn(length(t_noise), 1), 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, V_noise, 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