diff --git a/matlab/frf_data.mat b/matlab/frf_data.mat new file mode 100644 index 0000000..d04f1d5 Binary files /dev/null and b/matlab/frf_data.mat differ diff --git a/matlab/identif_analyze_all.m b/matlab/identif_analyze_all.m new file mode 100644 index 0000000..51f7159 --- /dev/null +++ b/matlab/identif_analyze_all.m @@ -0,0 +1,94 @@ +%% 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 +meas_data = {}; +for i = 1:6 + meas_data(i) = {load(sprintf('mat/frf_data_exc_strut_%i_noise_lf.mat', i), 't', 'Va', 'Vs', 'de')}; +% load(sprintf('mat/frf_data_exc_strut_%i_noise_hf.mat', i), 't', 'Va', 'Vs', 'de'); +end + +%% +Ts = (meas_data{1}.t(end) - (meas_data{1}.t(1)))/(length(meas_data{1}.t)-1); +Fs = 1/Ts; + +win = hanning(ceil(1*Fs)); % Hannning Windows + +%% DVF +[~, f] = tfestimate(meas_data{1}.Va, meas_data{1}.de, win, [], [], 1/Ts); + +G_dvf = zeros(length(f), 6, 6); +for i = 1:6 + G_dvf(:, :, i) = tfestimate(meas_data{i}.Va, meas_data{i}.de, win, [], [], 1/Ts); +end + +%% IFF +[~, f] = tfestimate(meas_data{1}.Va, meas_data{1}.de, win, [], [], 1/Ts); + +G_iff = zeros(length(f), 6, 6); +for i = 1:6 + G_iff(:, :, i) = tfestimate(meas_data{i}.Va, meas_data{i}.Vs, win, [], [], 1/Ts); +end + +%% +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i =1:6 + plot(f, abs(G_dvf(:,i, 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; +hold on; +for i =1:6 + plot(f, 180/pi*angle(G_dvf(:,i, 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([5, 5e3]); + +%% +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i =1:6 + plot(f, abs(G_iff(:,i, 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; +hold on; +for i =1:6 + plot(f, 180/pi*angle(G_iff(:,i, 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([5, 5e3]); \ No newline at end of file diff --git a/matlab/identif_save.m b/matlab/identif_save.m index a228022..6fe9afa 100644 --- a/matlab/identif_save.m +++ b/matlab/identif_save.m @@ -24,7 +24,7 @@ t = data(:, end); % Time [s] % And we save this to a =mat= file: -strut_number = 1; +strut_number = 6; % save(sprintf('mat/frf_data_exc_strut_%i_noise.mat', strut_number), 't', 'Va', 'Vs', 'de'); % save(sprintf('mat/frf_data_exc_strut_%i_noise_lf.mat', strut_number), 't', 'Va', 'Vs', 'de'); diff --git a/matlab/src/generateShapedNoise.m b/matlab/src/generateShapedNoise.m new file mode 100644 index 0000000..51748c2 --- /dev/null +++ b/matlab/src/generateShapedNoise.m @@ -0,0 +1,53 @@ +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 - [-] + +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 + +t_noise = 0:args.Ts:args.exc_duration; + +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 + +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, V_noise, V_smooth_end]; +t_exc = args.Ts*[0:1:length(V_exc)-1]; + +U_exc = [t_exc; 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/matlab/src/generateSweepExc.m b/matlab/src/generateSweepExc.m new file mode 100644 index 0000000..252b8c5 --- /dev/null +++ b/matlab/src/generateSweepExc.m @@ -0,0 +1,76 @@ +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 - [-] + +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 + +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 + +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 + +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 + +V_exc = [V_smooth_start, V_sweep, V_smooth_end]; +t_exc = args.Ts*[0:1:length(V_exc)-1]; + +U_exc = [t_exc; V_exc];