% Série de Fourier de uma onda quadrada

clc
clearvars
close all

% tolerância
tol=1e-12;

% período e frequência da função
T=5;
om1=2*pi/T;

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

% 1o. ponto - série de Fourier da função (forma real)
% vetor de tempo
np=500;
t=linspace(-T/2,T/2,np);
% função original
go=zeros(1,length(t));
for it=1:length(t)
  if t(it) < 0
    go(it)=-1;
  else
    go(it)=1;
  end
end
% número de termos oscilatórios
nto=5; % 1, 2, 3, 5, 11, 25, 51
% série de Fourier (forma real)
if nto > 0
  for j=1:nto
    om(j)=j*om1; 
    A(j)=0;
    B(j)=(2/(j*pi))*(1-cos(j*pi));
    G(j)=sqrt(A(j)^2+B(j)^2);
    if G(j) < tol
        fi(j)=0;
    else
        fi(j)=atan2(A(j),B(j));
    end
    gsr1(j,:)=A(j)*cos(j*om1*t)+B(j)*sin(j*om1*t); % forma 1
    gsr2(j,:)=G(j)*sin(j*om1*t+fi(j)); % forma 2
  end
end
% visualização da função e de sua série de Fourier
figure(1)
plot(t,go,'b.',t,sum(gsr1,1),'r-')
hold on
plot([-T/2 T/2],[0 0],'k-',[0 0],[-2 2],'k-')
hold off
axis([-T/2 T/2 -2 2])
xlabel('t(s)')
ylabel('g(t)')
grid on

% % 2o. ponto - série de Fourier da função (forma complexa)
% % montagem da série de Fourier 
% Gb=[]; gb=[]; psib=[];
% if nto > 0
%   for j=-nto:1:nto
%     % coeficientes
%     if j == 0
%       Gbaux=0;
%     else
%       Gbaux=(-1i/(j*pi))*(1-cos(j*pi));
%     end
%     % módulos de coeficientes
%     gaux=abs(Gbaux);
%     % argumentos de coeficientes
%     if gaux < tol
%       psiaux=0;
%     else
%       psiaux=angle(Gbaux);
%     end
%     % montagem de vetores de coeficientes, módulos e argumentos
%     Gb=[Gb; Gbaux];
%     gb=[gb; gaux];
%     psib=[psib; psiaux];
%   end
%   omb=[-om(length(om):-1:1) 0 om];
%   % série de Fourier (forma complexa)
%   for j=1:length(omb)
%     gsc(j,:)=Gb(j)*exp(i*omb(j)*t);
%   end
% end
% % visualização da função e de sua série de Fourier
% % em formas real e complexa
% figure(2)
% plot(t,go,'b.',t,sum(gsr1,1),'r-')
% hold on
% plot([-T/2 T/2],[0 0],'k-',[0 0],[-2 2],'k-')
% hold off
% axis([-T/2 T/2 -2 2])
% xlabel('t(s)')
% ylabel('g(t)')
% grid on
% hold on
% pause(3)
% plot(t,sum(gsc,1),'g-')
% pause(3)
% hold off

% % 3o. ponto - representação da função nas perspectivas
% % do tempo e da frequência (5 componentes)
% % (antes, atribuir nto = 5, na linha 33)
% umv=ones(size(t));
% figure(3)
% plot3(t,om1*umv,gsr2(1,:),'c',t,2*om1*umv,gsr2(2,:),'g',...
%       t,3*om1*umv,gsr2(3,:),'r',t,4*om1*umv,gsr2(4,:),'b',...
%       t,5*om1*umv,gsr2(5,:),'m',...
%       t,-om1*umv,sum(gsr2,1),'k');
% grid on  
% xlabel('t(s)');
% ylabel('\omega (rad/s)');
% zlabel('valor do sinal');
% %view(0,0); % visão no tempo
% %view(90,0); % visão na frequência

% % 4o. ponto - espectros unilaterais de amplitude e fase da função
% % (antes, atribuir a nto o número de termos desejado, tal que nto >= 1)
% figure(4)
% % espectro de amplitude
% subplot(2,1,1)
% stem(om,G,'b.')
% axis([0 ceil(max(om)) 0 1.5])
% xlabel('\omega (rad/s)');
% ylabel('amplitude');
% grid on
% % espectro de fase
% subplot(2,1,2)
% stem(om,fi,'b.')
% axis([0 ceil(max(om)) -pi +pi])
% xlabel('\omega (rad/s)');
% ylabel('fase (rad)');
% grid on
 
% % 5o. ponto - espectros bilaterais de amplitude e fase da função
% % via expressão analítica dos coeficientes complexos
% % visualização de espectros
% figure(5)
% % espectro de amplitude
% subplot(2,1,1)
% stem(omb,gb,'r.')
% hold on
% plot([0 0],[0 1],'k-')
% hold off
% axis([min(omb)-2 max(omb)+2 0 1]);
% xlabel('\omega (rad/s)');
% ylabel('amplitude');
% grid on
% % espectro de fase
% subplot(2,1,2)
% stem(omb,psib,'r.')
% hold on
% plot([0 0],[-pi pi],'k-')
% hold off
% axis([min(omb)-2 max(omb)+2 -pi +pi]);
% xlabel('\omega (rad/s)');
% ylabel('fase (rad)');
% grid on

% % 6o. ponto - espectros bilaterais de amplitude e fase da função
% % via relações com os espectros unilaterais correspondentes
% % determinação de componentes (amplitudes e fases)
% g2=G/2;
% g2b=[g2(length(g2):-1:1) 0 g2]; 
% for j=1:nto
%   if G(j) < tol
%     psi2(j)=0;
%   else
%     psi2(j)=fi(j)-(pi/2);
%   end
% end
% psi2b=[-psi2(length(psi2):-1:1) 0 psi2];
% % visualização de espectros
% figure(6)
% % espectro de amplitude
% subplot(2,1,1)
% stem(omb,gb,'r.')
% hold on
% plot([0 0],[0 1],'k-')
% hold off
% axis([min(omb)-2 max(omb)+2 0 1]);
% xlabel('\omega (rad/s)');
% ylabel('amplitude');
% grid on
% hold on 
% pause(2)
% stem(omb,g2b,'g.')
% pause(2)
% hold off
% % espectro de fase
% subplot(2,1,2)
% stem(omb,psib,'r.')
% hold on
% plot([0 0],[-pi pi],'k-')
% hold off
% axis([min(omb)-2 max(omb)+2 -pi +pi]);
% xlabel('\omega (rad/s)');
% ylabel('fase (rad)');
% grid on
% hold on 
% pause(2)
% stem(omb,psi2b,'g.')
% pause(2)
% hold off