clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=0.2;
ws=2*pi/T;
disp(['T = ',num2str(T)])
disp(['ws = ',num2str(ws)])
W=2*pi*6/(2*pi);
%%%%%%%% %%%%%%%%%%% tc=0:0.0001:10; %%% simulating the continous domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['W (max freq of the signal x(t)) = ',num2str(W)])
if ws/2>W 
    YesNo='Yes';
else
     YesNo='No';
end
disp(['Nyquist ws/2>W ? ',YesNo])
%xc=sin(W*tc);
%xs=sin(W*ts)+j*cos(W*ts);
ts=0:T:10;
xs=sin(W*ts)+0.5*randn(1,length(ts));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot(tc,xc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%hold on
figure
stem(ts,xs,'r')
title('DATA')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L=length(xs); %%% length of the signal
N=L; %%% you can choose N>=L 
FXs=fft(xs,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
stem((0:N-1),abs(FXs).^2)
hold on
plot([N/2 N/2],[0 max(abs(FXs).^2)],'r--','LineWidth',2)
set(gca,'fontSize',20,'fontWeight','Bold')
xlabel('k')
xlim([0 N-1])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
Omega=(0:N-1)*(2*pi/N);
stem(Omega,abs(FXs).^2)
hold on
plot([pi pi],[0 max(abs(FXs).^2)],'r--','LineWidth',2)
set(gca,'fontSize',20,'fontWeight','Bold')
xlabel('\Omega')
xlim([0 2*pi])

figure
w=0:(ws/2)/(N/2):ws-(ws/2)/(N/2);
stem(w,abs(FXs).^2)
hold on
plot([pi/T pi/T],[0 max(abs(FXs).^2)],'r--','LineWidth',2)
set(gca,'fontSize',20,'fontWeight','Bold')
xlabel('w=\Omega/T')
xlim([0 2*pi/T])
title(['T = ', num2str(T)])


figure
FXs2=FXs(1:fix(N/2));
w=0:(ws/2)/(N/2):ws/2-(ws/2)/(N/2);
stem(w,abs(FXs2),'r')
hold on
set(gca,'fontSize',20,'fontWeight','Bold')
xlabel('w=\Omega/T')
xlim([0 pi/T])
title(['T = ', num2str(T)])


xlim([0 ws/2])