% Análise de Fourier de uma onda quadrada periódica
% pelas vias numérica e analítica

clc
clearvars
close all

% tolerância
tol=1e-12;

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

% % 1o. ponto - montagem e visualização das funções
%   % função periódica original (um período)
% npt=1024;
% t=linspace(0,T,npt);
% go=zeros(1,length(t));
%   % função original
% for it=1:length(t)
%   if t(it) < (T/2)
%     go(it)=1;
%   else
%     go(it)=-1;
%   end
% end
%   % função periódica "amostrada" (um período)
%     % número de pontos amostrais
% N=32; % 32, 64 e 128 (em cômputo posterior da fase)
%     % intervalo amostral 
% dt=T/N;
%     % vetor de tempo correspondente
% t2=0:dt:(N-1)*dt;
%     % localização da descontinuidade
% Nd=floor(T/(2*dt));
%     % montagem da função
% gd1=ones(1,Nd);
% gd2=-ones(1,(N-Nd));
% gd=[gd1 gd2];
%   % visualização das funções original e "amostrada"
% figure(1)
% plot(t,go,'b.')
% axis([0 T -1.5 1.5])
% xlabel('t(s)');
% ylabel('g(t)');
% legend('original')
% grid on
% pause(3)
% hold on
% plot(t2,gd,'ro')
% legend('original','amostrada')
% plot([0 T],[0 0],'k-')
% hold off

% % 2o. ponto - cômputo da tdf e tidf, tb. via fft e ifft
% % (comentar linhas 41 a 54)
%   % cálculo da tdf
% Xb=tdf(gd);
% % Xb=fft(gd);
%   % cálculo da itdf
% gr=itdf(Xb);
% % gr=ifft(Xb);
%   % visualização
% figure(2)
% plot(t,go,'b.')
% axis([0 T -1.5 1.5])
% xlabel('t(s)');
% ylabel('g(t)');
% legend('original')
% grid on
% pause(3)
% hold on
% plot(t2,gd,'ro')
% legend('original','amostrada')
% pause(3)
% plot(t2,gr,'go')
% legend('original','amostrada','recuperada')
% plot([0 T],[0 0],'k-')
% hold off

% % 3o. ponto - espectros bilaterais de amplitude e fase da função
% % (comentar linhas 64 a 80)
%   % coeficientes complexos via expressão analítica 
% om=zeros(1,N);
% ga=[]; psia=[];
% for j=(-N/2+1):1:(N/2)
%   ij=j+(N/2);
%   om(ij)=j*om1;
%   if j == 0
%     Gbaux=0;
%   else
%     Gbaux= (-1i/(j*pi))*(1-cos(j*pi));
%   end
%   gaux=abs(Gbaux);
%   if gaux < tol
%     psiaux=0;
%   else
%     % psiaux=angle(Gbaux);
%     psiaux=atan2(imag(Gbaux),real(Gbaux));
%   end
%   ga=[ga gaux];
%   psia=[psia psiaux];
% end
%   % coeficientes complexos via tdf
% Gbn=Xb/N; % relação entre coeficientes e TDF
% gn=abs(Gbn);
% psin=atan2(imag(Gbn),real(Gbn));
%     % deslocamento para montagem de espectros bilaterais 
%     % em vista da periodicidade dos coeficientes numéricos
% gn=[gn((N/2+2):N) gn(1:N/2+1)];
% psin=[psin((N/2+2):N) psin(1:N/2+1)];
%   % visualização
% figure(3)
% g=[ga;gn]';
% h=stem(om,g);
% set(h(1),'Color','blue','Marker','o')
% set(h(2),'Color','red','Marker','s')
% legend('analítico','numérico')
% xlabel('\omega (rad/s)');
% ylabel('amplitude');
% grid on
% pause(3)
% figure(4)
% psi=[psia;psin]';
% h=stem(om,psi);
% set(h(1),'Color','blue','Marker','o')
% set(h(2),'Color','red','Marker','s')
% legend('analítico','numérico')
% xlabel('\omega (rad/s)');
% ylabel('fase (rad)');
% grid on