% Série de Fourier de uma função retangular com período crescente

clc
clearvars
close all

% dados da força
L=2; % 2s (passos 1 a 5); 1s e 0.1s (passo 6)
R=2/L; 
% período da força
% valor inicial: 4s
% valor final (máximo): 32s
T=4;

% frequência fundamental da força
om1=2*pi/T;

% número de termos oscilatórios
% valor inicial: 5 (passos 1 a 4)
% valor subsequente: 25 (início do passo 5)
% valor final (máximo): 200 (final do passo 5)
nto=5; 

% visualização de resultados preliminares
disp(' ')
disp(['A frequência fundamental da excitação é ',num2str(om1),' rad/s.'])
disp(' ')

% Após o 4o. passo, comentar linhas 48, 71, 72, 94, 109, 129, 130 e 138

% 1o. passo - série de Fourier, forma real, da força
% definição do vetor de tempo
tl=32;
np=40*tl;
t=linspace(-tl,tl,np); 
% composição e visualização da série
  % termo constante
A02=R*L/T;
F0=A02;
f0=F0*ones(size(t));
ftr=f0;
figure(1)
plot(t,ftr,'-b')
axis([-40 40 -0.5 R+0.5])
grid on
xlabel('t (s)')
ylabel('f(t) (N)')
pause(4) % pausa
  % termos oscilatórios
if nto > 0
  for j=1:nto
    A(j)=(2*R*L/T)*(sin(j*pi*L/T)/(j*pi*L/T));
    B(j)=0;
    F(j)=abs(A(j));
    if F(j) < 1e-10
      fif(j)=0;
    else
      if A(j) > 0
        fif(j)=pi/2;
      else
        fif(j)=-pi/2;
      end
    end
    om(j)=j*om1;
    ftr=ftr+F(j)*sin(om(j)*t+fif(j));
    plot(t,ftr,'-b')
    axis([-40 40 -0.5 R+0.5])
    grid on
    xlabel('t (s)')
    ylabel('f(t) (N)')
    legend(['nto = ',num2str(j)])
    pause(3) % pausa
  end
end  

% 2o. passo - espectros unilaterais de amplitude e fase da força
% (atribuir a nto, na linha 22, o número de termos desejado, com nto > 0)
figure(2)
subplot(2,1,1)
FU=[F0 F];
omu=[0 om];
stem(omu,FU,'.b') %. %:.
axis([0 1.1*max(omu) 0 1.1*max(FU)])
grid on
xlabel('\omega (rad/s)');
ylabel('amplitude (N)');
subplot(2,1,2)
fifu=[0 fif];
stem(omu,fifu,'.b') %. %:.
axis([0 1.1*max(omu) -pi +pi])
grid on
xlabel('\omega (rad/s)');
ylabel('fase (rad)');
pause(3) % pausa

% 3o. passo - série de Fourier da força, forma complexa, via relações
%             com coeficientes da forma real
% composição e visualização da série
  % termo constante
F0=R*L/T;
f0=F0*ones(size(t));
ftc=f0;
figure(3)
plot(t,ftc,'-r')
axis([-40 40 -0.5 R+0.5])
grid on
xlabel('t (s)')
ylabel('f(t) (N)')
pause(3) % pausa
  % termos oscilatórios
if nto > 0
  for j=1:nto
    Fb(j)=(R*L/T)*(sin(j*pi*L/T)/(j*pi*L/T));
    f(j)=abs(Fb(j));
    if f(j) < 1e-10
      psif(j)=0;
    else
      psif(j)=fif(j)-pi/2;
    end
    om(j)=j*om1;
    ftc=ftc+...
        f(j)*exp(1i*(-psif(j)))*exp(1i*(-om(j))*t)+...
        f(j)*exp(1i*psif(j))*exp(1i*om(j)*t);
    plot(t,ftc,'-r')
    axis([-40 40 -0.5 R+0.5])
    grid on
    xlabel('t (s)')
    ylabel('f(t) (N)')
    legend(['nto = ',num2str(j)])
    pause(3) % pausa
  end
end
hold on
  % superposição da série anterior (real)
plot(t,ftr,'b-')
legend('forma complexa','forma real')
hold off
pause(3) % pausa

% 4o. passo - espectros bilaterais de amplitude e fase da força
% (atribuir a nto, na linha 22, o número de termos desejado, com nto > 0)
figure(4)
subplot(2,1,1)
FB=[f F0 f];
omb=[-om 0 om];
stem(omb,FB,'.r') %. %:.
axis([1.1*min(omb) 1.1*max(omb) 0 1.1*max(FB)])
grid on
xlabel('\omega (rad/s)');
ylabel('Amplitude (N)');
subplot(2,1,2)
psib=[-psif 0 psif];
stem(omb,psib,'.r') %. %:.
axis([1.1*min(omb) 1.1*max(omb) -1.1*pi 1.1*pi])
grid on
xlabel('\omega (rad/s)');
ylabel('Fase (rad)');
 
% 5o. passo - quando T tende para infinito
% -> comentar linhas 48, 71, 72, 94, 109, 129, 130 e 138
% -> rodar para os seguintes valores de T e nto: 4/25, 8/50, 16/100 e 32/200
%    (dobrar nto ao dobrar T, para observar a mesma faixa de frequência)

% 6o. passo - influência da duração do pulso
% -> L em 1s e 0.1s, além de T/nto em 32/200
%    (produto RL constante, para manter valor inicial no espectro)