1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
|
%% script de test harmonique
clc;
% clear all;
F=reshape(repmat(cumprod([1e-2 repmat(10,1,5)]),9,1),1,9*6) .* repmat(1:9,1,6) ;
% F contient les fréquence 1e-2 2e-2 .. 1e-1 2e-1 .... 9e3
AllTest=cell(length(F),1) ;
simopt_begin = simget('test_ordre1');
open_system('test_ordre1');
set_param('test_ordre1/SinusInput','Frequency','frequency');
set_param('test_ordre1/SinusInput','Amplitude','1');
set_param('test_ordre1/SinusInput','Bias','0');
set_param('test_ordre1/SinusInput','Phase','0');
set_param('test_ordre1/SinusInput','Offset','0');
set_param('test_ordre1/SinusInput','SampleTime','0');
set_param('test_ordre1/SinusInput','SineType','Time based');
set_param('test_ordre1/SinusInput','TimeSource','Use simulation time');
disp('TEST HARMONIQUE D''IDENTIFICATION DE FONCTION DE TRANSFERT (SIMULATION)');
for ind_freq = 1:length(F)
OneTest.freq = F(ind_freq) ;
frequency = F(ind_freq) ;
simopt=simset(simget('test_ordre1'),'MaxStep',pi/(250*frequency));
sim('test_ordre1',[0 4*pi/frequency],simopt);
OneTest.Scope = ScopeDataFrequentielle ;
AllTest{ind_freq} = OneTest ;
disp(sprintf('Test frequence %4.4f Hz OK',frequency));
end;
sim('test_ordre1',[0 0],simopt_begin);
disp('TEST HARMONIQUE D''IDENTIFICATION DE FONCTION DE TRANSFERT (ANALYSE)');
Gain = zeros(1,length(AllTest)) ;
Dephase = zeros(1,length(AllTest)) ;
Frequence = F ;
for ind_freq = 1:length(AllTest)
frequency = AllTest{ind_freq}.freq ;
t = AllTest{ind_freq}.Scope.time ;
y_in = AllTest{ind_freq}.Scope.signals.values(:,1) ;
y_out = AllTest{ind_freq}.Scope.signals.values(:,2) ;
[ymin_in,indmin_in] = min(y_in) ;
[ymax_in,indmax_in] = max(y_in(1:indmin_in)) ;
[ymax2_in,indmax2_in] = max(y_in(indmin_in:end)) ;
indmax2_in = indmax2_in+indmin_in-1;
[ymin_out,indmin_out] = min(y_out(indmin_in:end)) ;
indmin_out = indmin_out + indmin_in - 1;
[ymax_out,indmax_out] = max(y_out(indmax_in:indmax2_in)) ;
indmax_out = indmax_out + indmax_in - 1 ;
Gain(ind_freq) = 0.5*(ymax_out-ymin_out)/(ymax_in-ymin_in);
Dephase(ind_freq) = (t(indmax_out)-t(indmax_in))*frequency ;
end;
figure;
subplot(2,1,1)
semilogx(Frequence,20*log10(abs(Gain/2)));
title(['Amplitude H(jw)']);
ylabel('Amplitude [dB]');
xlabel('W [rad/s]');
grid on;
subplot(2,1,2);
semilogx(Frequence,Dephase*180/pi);
title('Phase H(jw)');
ylabel('Angle [deg.]');
xlabel('W [rad /s]');
grid on; |