Analyze first measurements
This commit is contained in:
		@@ -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]);
 | 
			
		||||
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
 | 
			
		||||
 
 | 
			
		||||
@@ -1,4 +1,4 @@
 | 
			
		||||
% Save Data
 | 
			
		||||
% =frf_save.m= - Save Data
 | 
			
		||||
% :PROPERTIES:
 | 
			
		||||
% :header-args: :tangle matlab/frf_save.m
 | 
			
		||||
% :END:
 | 
			
		||||
 
 | 
			
		||||
@@ -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');
 | 
			
		||||
 
 | 
			
		||||
							
								
								
									
										54
									
								
								matlab/src/generateSinIncreasingAmpl.m
									
									
									
									
									
										Normal file
									
								
							
							
						
						
									
										54
									
								
								matlab/src/generateSinIncreasingAmpl.m
									
									
									
									
									
										Normal file
									
								
							@@ -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];
 | 
			
		||||
		Reference in New Issue
	
	Block a user