measure LF FRF

matlab/identif_analyze.m

@@ -9,7 +9,9 @@ addpath('./src/');
 % Test with one APA
 
 %% Load measurement data for APA number 1
-load(sprintf('mat/frf_data_%i_sweep_lf.mat', 2), 't', 'Va', 'Vs', 'de', 'da');
+strut_number = 1;
+% load(sprintf('mat/frf_data_exc_strut_%i_noise_lf.mat', strut_number), 't', 'Va', 'Vs', 'de');
+load(sprintf('mat/frf_data_exc_strut_%i_noise_hf.mat', strut_number), 't', 'Va', 'Vs', 'de');
 
 % Compute transfer functions:
 
@@ -18,31 +20,17 @@ Fs = 1/Ts;
 
 win = hanning(ceil(1*Fs)); % Hannning Windows
 
+%% DVF
 [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;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
-plot(f, coh_dvf);
-plot(f, coh_d);
-plot(f, coh_iff);
-hold off;
-set(gca, 'XScale', 'log');
 
-%%
+ax1 = nexttile([2,1]);
-figure;
-tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
-
-ax1 = nexttile;
 hold on;
-plot(f, abs(G_dvf));
+for i =1:6
-plot(f, abs(G_d));
+    plot(f, abs(G_dvf(:,i)));
+end
 hold off;
 set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 ylabel('Amplitude $V_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]);
@@ -50,8 +38,9 @@ hold off;
 
 ax2 = nexttile;
 hold on;
-plot(f, 180/pi*angle(G_dvf));
+for i =1:6
-plot(f, 180/pi*angle(G_d));
+    plot(f, 180/pi*angle(G_dvf(:,i)));
+end
 hold off;
 set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
 xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
@@ -61,29 +50,58 @@ yticks(-360:90:360);
 linkaxes([ax1,ax2],'x');
 xlim([5, 5e3]);
 
-figure;
+%% IFF
-tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+[G_iff, f] = tfestimate(Va, Vs, win, [], [], 1/Ts);
 
-ax1 = nexttile;
+figure;
-plot(f, abs(G_iff));
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+for i =1:6
+    plot(f, abs(G_iff(:,i)));
+end
+hold off;
 set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 ylabel('Amplitude $V_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]);
 hold off;
 
 ax2 = nexttile;
-plot(f, 180/pi*angle(G_iff));
+hold on;
+for i =1:6
+    plot(f, 180/pi*angle(G_iff(:,i)));
+end
+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([0.1, 10]);
+xlim([5, 5e3]);
 
-% Comparison of all APA
+%%
+[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts);
+[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts);
 
-%% Load all the measurements
+%%
-meas_data = {};
+figure;
-for i = 1:7
+hold on;
-    meas_data(i) = {load(sprintf('mat/frf_data_%i.mat', i), 't', 'Va', 'Vs', 'de', 'da')};
+for i =1:6
+    plot(f, coh_dvf(:,i));
 end
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Coherence');
+ylim([0,1]);
+
+%%
+figure;
+hold on;
+for i =1:6
+    plot(f, coh_iff(:,i));
+end
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Coherence');
+ylim([0,1]);

matlab/identif_save.m

@@ -17,19 +17,17 @@ close(f);
 
 data = SimulinkRealTime.utils.getFileScopeData('data/data.dat').data;
 
-da = data(:, 1); % Excitation Voltage (input of PD200) [V]
+de = data(:, 1:6); % Measurment displacement (encoder) [m]
-de = data(:, 2); % Measured voltage (force sensor) [V]
+Vs = data(:, 7:12); % Measured voltage (force sensor) [V]
-Vs = data(:, 3); % Measurment displacement (encoder) [m]
+Va = data(:, 13); % Excitation Voltage [V]
-Va = data(:, 4); % Measurement displacement (attocube) [m]
 t  = data(:, end); % Time [s]
 
 
-
 % And we save this to a =mat= file:
-apa_number = 1;
+strut_number = 1;
-% leg_number = 4;
 
-save(sprintf('mat/frf_data_leg_coder_%i_noise.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da');
+% save(sprintf('mat/frf_data_exc_strut_%i_noise.mat', strut_number), 't', 'Va', 'Vs', 'de');
-% save(sprintf('mat/frf_data_leg_coder_%i_sweep.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da');
+% save(sprintf('mat/frf_data_exc_strut_%i_noise_lf.mat', strut_number), 't', 'Va', 'Vs', 'de');
-% save(sprintf('mat/frf_data_leg_coder_%i_noise_hf.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da');
+% save(sprintf('mat/frf_data_exc_strut_%i_sweep.mat', strut_number), 't', 'Va', 'Vs', 'de');
-% save(sprintf('mat/frf_data_leg_coder_%i_add_mass_closed_circuit.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da');
+save(sprintf('mat/frf_data_exc_strut_%i_noise_hf.mat', strut_number), 't', 'Va', 'Vs', 'de');
+% save(sprintf('mat/frf_data_exc_strut_%i_add_mass_closed_circuit.mat', strut_number), 't', 'Va', 'Vs', 'de');

matlab/identif_setup.m

@@ -7,8 +7,8 @@ s = zpk('s');
 addpath('./src/');
 
 %% Simulation configuration
-Fs   = 10e3; % Sampling Frequency [Hz]
+Fs = 10e3; % Sampling Frequency [Hz]
-Ts   = 1/Fs; % Sampling Time [s]
+Ts = 1/Fs; % Sampling Time [s]
 
 %% Data record configuration
 Trec_start = 5;  % Start time for Recording [s]
@@ -27,7 +27,7 @@ V_noise = generateShapedNoise('Ts', 1/Fs, ...
 %% Sweep Sine
 gc = 0.1;
 xi = 0.5;
-wn = 2*pi*94.3;
+wn = 2*pi*92.7;
 
 % 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);
@@ -41,8 +41,17 @@ V_sweep = generateSweepExc('Ts',      Ts, ...
                            'sweep_type',   'log', ...
                            'V_exc', G_sweep*1/(1 + s/2/pi/500));
 
+V_sweep_lf = generateSweepExc('Ts',      Ts, ...
+                           'f_start', 0.1, ...
+                           'f_end',   10, ...
+                           'V_mean',  3.25, ...
+                           't_start', Trec_start, ...
+                           'exc_duration', Trec_dur, ...
+                           'sweep_type',   'log', ...
+                           'V_exc', 0.2);
+
 %% High Frequency Shaped Noise
-[b,a] = cheby1(10, 2, 2*pi*[300 2e3], 'bandpass', 's');
+[b,a] = cheby1(10, 2, 2*pi*[240 2e3], 'bandpass', 's');
 wL = 0.005*tf(b, a);
 
 V_noise_hf = generateShapedNoise('Ts', 1/Fs, ...
@@ -52,6 +61,17 @@ V_noise_hf = generateShapedNoise('Ts', 1/Fs, ...
                               'smooth_ends', true, ...
                               'V_exc', wL);
 
+%% Low Frequency Shaped Noise
+[b,a] = cheby1(10, 2, 2*pi*[10 260], 'bandpass', 's');
+wL = 0.005*tf(b, a);
+
+V_noise_lf = generateShapedNoise('Ts', 1/Fs, ...
+                              'V_mean', 3.25, ...
+                              't_start', Trec_start, ...
+                              'exc_duration', Trec_dur, ...
+                              'smooth_ends', true, ...
+                              'V_exc', wL);
+
 %% Sinus excitation with increasing amplitude
 V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ...
                                   'V_mean', 3.25, ...
@@ -61,20 +81,33 @@ V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ...
                                   't_start', Trec_start, ...
                                   'smooth_ends', true);
 
+%% Zero Excitation
+% Trec_start = 10;  % Start time for Recording [s]
+% Trec_dur   = 10; % Recording Duration [s]
+% 
+% Tsim = 2*Trec_start + Trec_dur;  % Simulation Time [s]
+
+V_zero = generateShapedNoise('Ts', 1/Fs, ...
+                             'V_mean', 3.25, ...
+                             't_start', Trec_start, ...
+                             'exc_duration', Trec_dur, ...
+                             'smooth_ends', true, ...
+                             'V_exc', tf(0));
+
 %% Select the excitation signal
-V_exc = timeseries(V_noise(2,:), V_noise(1,:));
+V_exc = timeseries(V_noise_hf(2,:), V_noise_hf(1,:));
 
 %% Plot
 figure;
 tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');
 
 ax1 = nexttile;
-plot(V_exc(1,:), V_exc(2,:));
+plot(V_exc.Time, squeeze(V_exc.Data));
 xlabel('Time [s]'); ylabel('Amplitude [V]');
 
 ax2 = nexttile;
-win = hanning(floor(length(V_exc)/8));
+win = hanning(floor(length(squeeze(V_exc.Data))/8));
-[pxx, f] = pwelch(V_exc(2,:), win, 0, [], Fs);
+[pxx, f] = pwelch(squeeze(V_exc.Data), win, 0, [], Fs);
 plot(f, pxx)
 xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]');
 set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');