% FFT de vibração harmônica

clc
clearvars
close all

% tolerância
tol=1e-10;

% criação de vibração harmônica
f1=20; % frequência do sinal (Hz)
fs=100; % frequência de amostragem: fs >= 2.56*f1 (Hz)
T=10; % tempo de aquisição (s)
t=0:1/fs:T-1/fs; % vetor de tempo (s)
x1=sin(2*pi*f1*t); % vibração (10^6m)

% cômputo da transformada de Fourier
N=length(t); % comprimento do vetor de tempo
X1=fft(x1); % aplicação da FFT
f=(1/T)*(0:N-1); % vetor de frequência (Hz)

% visualização de amplitude e fase (f > 0)
figure(1);
subplot(1,2,1)
plot(f,abs(X1)/N)
grid on
xlabel('frequência (Hz)'); ylabel('amplitude (\mum.s)')
title('Espectro de amplitude')
axis([0 fs/2 0 0.55]) % limitação até fs/2
subplot(1,2,2)
if real(X1) < tol
  inz=find(abs(imag(X1)) > tol);
  X1c=zeros(size(X1))+1i*zeros(size(X1));
  X1c(inz)=X1(inz);
  plot(f,atan2(imag(X1c),real(X1c)))
elseif imag(X1) < tol
  inz=find(abs(real(X1)) > tol);
  X1c=zeros(size(X1))+1i*zeros(size(X1));
  X1c(inz)=X1(inz);
  plot(f,atan2(imag(X1c),real(X1c)))
else
  plot(f,atan2(imag(X1),real(X1)))
end
grid on
xlabel('frequência (Hz)'); ylabel('fase (rad)')
title('Espectro de fase')
axis([0 fs/2 -pi pi]) % limitação até fs/2