function [t,x,xh,xe,w,f,X] = ajusta(arquivo,t1,t2,t3,t4,t5,t6,imprime)
% [t,x,xh,xe,w,f,X] = ajusta(arquivo,t1,t2,t3,t4,t5,t6,imprime)
%
% Análise de sinal de vibração amortecida
%
% t = vetor de tempo [n,1];
% x = vetor de velocidade angular [n,1];
% xh = reta média ajustada à envoltória do sinal no intervalo [t1:t2];
% xe = exponencial média ajustada à envoltória do sinal no intervalo [t3:t4];
% w = frequencia instantânea ajustada no intervalo [t5:t6];
% f = vetor de frequencia [n/2,1];
% X = espectro do sinal [n/2,1];
% arquivo = nome do arquino no formato csv
% t1 = início do intervalo para ajuste da reta à envoltória;
% t2 = termino do intervalo para ajuste da reta à envoltória;
% t3 = início do intervalo para ajuste da exponencial à envoltória;
% t4 = termino do intervalo para ajuste da exponencial à envoltória;
% t5 = início do intervalo para ajuste da reta à fase instantânea;
% t6 = termino do intervalo para ajuste da reta à fase instantênea;
% imprime = se existir as figuras são plotadas.
%
% Walter Ponge-Ferreira
% Cotia, 02.02.2021
%


if nargin<1,
  arquivo = 'giroscopio01';
endif
if nargin<2,
  t1 = 0.3;
endif
if nargin<3,
  t2 = 9;
endif
if nargin<4,
  t3 = 8;
endif
if nargin<5,
  t4 = 14.5;
endif
if nargin<6,
  t5 = 0.0;
endif

% Leitura do sinal
pkg load signal

eval(sprintf('load %s.csv -ascii',arquivo))
eval(sprintf('t = %s(:,1);',arquivo))
eval(sprintf('x = %s(:,2);',arquivo))
eval(sprintf('y = %s(:,3);',arquivo))
eval(sprintf('z = %s(:,4);',arquivo))
eval(sprintf('w = %s(:,5);',arquivo))
if nargin<7,
  t6 = length(t)*diff(t(1:2));
endif

dt = diff(t(1:2));
n = 2^floor(log2(length(t)));
t = t(1:n);
x = x(1:n);
y = y(1:n);
z = z(1:n);
w = w(1:n);

% Registro do sinal
figure(1)
clf
plot(t,x)
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
if nargin>=8,
  print -deps fig01.eps
  print -djpeg fig01.jpg
endif

% Determinação do espectro
X = fft(x)/n*2;
T = n*dt;
df = 1/T;
f = [0:n/2-1]'*df;
[Xmax,i] = max(abs(X));
figure(2)
clf
plot(f(1:128),abs(X(1:128)),f(i),abs(X(i)),'*r')
text(f(i),abs(X(i)),sprintf('%5.2f Hz',f(i)))
xlabel('frequencia - f / Hz')
ylabel('velocidade rotacao - |\Omega| / rad/s')
grid
if nargin>=8,
  print -deps fig02.eps
  print -djpeg fig02.jpg
endif

% Determinação da Envoltória
xx = hilbert(x);
j = find(t>t1&t<t2);
figure(3)
plot(t,abs(xx),t(j),abs(xx(j)),'r')
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
shg

% Ajusto de reta
P = polyfit(t(j),abs(xx(j)),1);
xh = polyval(P,t);
plot(t,abs(xx),t(j),abs(xx(j)),t,xh,'k')
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
axis([0,16,0,0.8])
shg
if nargin>=8,
  print -deps fig03.eps
  print -djpeg fig03.jpg
endif

figure(4)
clf
plot(t,x,t,xh,t,-xh)
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
shg
if nargin>=8,
  print -deps fig04.eps
  print -djpeg fig04.jpg
endif

% Ajuste de exponencial
j2 = find(t>t3&t<t4);
Pe = polyfit(t(j2),log(abs(xx(j2))),1);
xe = exp(Pe(2))*exp(Pe(1)*t);
figure(5)
clf
semilogy(t,abs(xx),t(j2),xe(j2),'r')
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
shg
if nargin>=8,
  print -deps fig05.eps
  print -djpeg fig05.jpg
endif

figure(6)
clf
plot(t,abs(xx),t(j),abs(xx(j)),t(j),xh(j),'k',t(j2),xe(j2),'r')
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
shg
if nargin>=8,
  print -deps fig06.eps
  print -djpeg fig06.jpg
endif

figure(7)
clf
plot(t,x,t(j),xh(j),'k',t(j),-xh(j),'k',t(j2),xe(j2),'r',t(j2),-xe(j2),'r')
xlabel('tempo - t / s')
ylabel('velocidade rotacao - \omega / rad/s')
grid
shg
if nargin>=8,
  print -deps fig07.eps
  print -djpeg fig07.jpg
endif

% Determinação da fase e frequência instantânea
o = angle(xx);
o = unwrap(o);
j3 = find(t>t5&t<t6);
Po = polyfit(t(j3),o(j3),1);
oo = polyval(Po,t);
figure(8)
clf
plot(t,o,t,oo,':r',t([j3(1),j3(end)]),oo([j3(1),j3(end)]),'*r')
xlabel('tempo - t / s')
ylabel('fase - \theta / rad')
text(0.6*T,0.6*max(o),sprintf('f_o = %5.2f Hz',Po(1)/2/pi))
grid
shg
if nargin>=8,
  print -deps fig08.eps
  print -djpeg fig08.jpg
endif

figure(9)
clf
ww = Po(1)*ones(n,1);
w = diff(o)./diff(t);
w = [w(1);w];
plot(t,w/2/pi,t,ww/2/pi,'r')
xlabel('tempo - t / s')
ylabel('frequência - f / Hz')
text(0.6*T,max(w/2/pi),sprintf('f_o = %5.2f Hz',Po(1)/2/pi))
grid
shg
if nargin>=8,
  print -deps fig09.eps
  print -djpeg fig09.jpg
endif

endfunction