diff --git a/matlab/frf_analyze.m b/matlab/frf_analyze.m index 15c30c6..c497aac 100644 --- a/matlab/frf_analyze.m +++ b/matlab/frf_analyze.m @@ -5,7 +5,89 @@ %% Load all the measurements -meas_data = {}; -for i = 1:7 - meas_data(i) = {load(sprintf('mat/frf_data_%i.mat', i), 't', 'Va', 'Vs', 'd')}; -end +% 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') + +%% +figure; +plot(t, de); + +%% +figure; +plot(t, Va); + +%% +Ts = (t(end) - t(1))/(length(t)-1); +Fs = 1/Ts; + +win = hanning(ceil(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'); + +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; + +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]); + +%% +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; + +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]); \ No newline at end of file diff --git a/matlab/frf_measure.slx b/matlab/frf_measure.slx new file mode 100644 index 0000000..7d1b7ec Binary files /dev/null and b/matlab/frf_measure.slx differ diff --git a/matlab/frf_save.m b/matlab/frf_save.m index 6b5830f..1bf41ec 100644 --- a/matlab/frf_save.m +++ b/matlab/frf_save.m @@ -17,10 +17,10 @@ close(f); 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] @@ -29,4 +29,4 @@ t = data(:, end); % Time [s] 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'); diff --git a/matlab/frf_setup.m b/matlab/frf_setup.m index 568cf4a..c5072a1 100644 --- a/matlab/frf_setup.m +++ b/matlab/frf_setup.m @@ -1,20 +1,30 @@ +s = tf('s'); +addpath('src') + +%% Fs = 10e3; % Sampling Frequency [Hz] Ts = 1/Fs; % Sampling Time [s] -Tsim = 110; % Simulation Time [s] - Trec_start = 5; % Start time for Recording [s] Trec_dur = 100; % Recording Duration [s] +Tsim = 2*Trec_start + Trec_dur; % Simulation Time [s] + %% Sweep Sine +gc = 0.1; +xi = 0.5; +wn = 2*pi*94.3; + +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, ... + 'f_end', 2e3, ... 'V_mean', 3.25, ... '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)); %% Shaped Noise V_noise = generateShapedNoise('Ts', 1/Fs, ... @@ -22,22 +32,25 @@ V_noise = generateShapedNoise('Ts', 1/Fs, ... 't_start', Trec_start, ... 'exc_duration', Trec_dur, ... 'smooth_ends', true, ... - 'V_exc', 0.05/(1 + s/2/pi/10)); + 'V_exc', 0.00/(1 + s/2/pi/50)); %% Select the excitation signal -V_exc = V_noise; +V_exc = timeseries(V_noise(2,:), V_noise(1,:)); 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)); -[pxx, f] = pwelch(V_exc(2,:), win, 0, [], Fs); +win = hanning(floor(length(V_exc.Data)/8)); +[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'); xlim([1, Fs/2]); ylim([1e-10, 1e0]); + +%% Save +save('./frf_data.mat', 'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', 'V_exc'); diff --git a/matlab/mat/frf_data_1_huddle.mat b/matlab/mat/frf_data_1_huddle.mat new file mode 100644 index 0000000..2f42741 Binary files /dev/null and b/matlab/mat/frf_data_1_huddle.mat differ diff --git a/matlab/mat/frf_data_1_noise.mat b/matlab/mat/frf_data_1_noise.mat new file mode 100644 index 0000000..716aa95 Binary files /dev/null and b/matlab/mat/frf_data_1_noise.mat differ diff --git a/matlab/mat/frf_data_1_sweep.mat b/matlab/mat/frf_data_1_sweep.mat new file mode 100644 index 0000000..433eb3c Binary files /dev/null and b/matlab/mat/frf_data_1_sweep.mat differ