clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=0.01;
ws=2*pi/T;
disp(['T = ',num2str(T)])
disp(['ws = ',num2str(ws)])
W1=20; %%% this is an angular frequency (rad/sec)
W2=150; %%% this is an angular frequency (rad/sec)
W=max([W1,W2]); %%% ancho de banda
%%%%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])
disp(['Max Freq that I can see = ',num2str(ws/2),' rad/sec'])
disp(['Max Freq that I can see = ',num2str(ws/(4*pi)),' Hz'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%xc=sin(W*tc);
%xs=sin(W*ts)+j*cos(W*ts);
ts=0:T:10;
xs=2*sin(W1*ts)+cos(W2*ts);
figure
plot(ts,xs,'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L=length(xs);
N=L; %%% you can choose N>=L
%%%%%%%%%%%%%%%%%
FXs=fft(xs,N);
%%%%%%%%%%%%%%%%
figure
w=(0:N-1)*(2*pi/(N*T));
stem(w,abs(FXs))
hold on
plot([pi/T pi/T],[0 max(abs(FXs))],'r--','LineWidth',2)
xlabel('w - (rad/sec)')
set(gca,'FontSize',20)
set(gca,'FontWeight','Bold')
%%%%%%%
figure
f=w/(2*pi); %%%% freq in Hz
plot(f,abs(FXs),'o-')
hold on
plot([pi/(T*2*pi) pi/(T*2*pi)],[0 max(abs(FXs))],'r--','LineWidth',2)
xlabel('f - (Hz)')
set(gca,'FontSize',20)
set(gca,'FontWeight','Bold')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ts1=0:T:5;
ts2=5+T:T:10;
xs2=2*sin(W1*ts1);
xs2(end+1:L)=cos(W2*ts2);
figure
plot(ts,xs2,'r')
xlabel('xs2')
set(gca,'FontSize',20)
set(gca,'FontWeight','Bold')
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%
FXs2=fft(xs2,N);
%%%%%%%%%%%%%%%%%%%%%%%
figure
f=w/(2*pi); %%%% freq in Hz
plot(f,abs(FXs2),'o-')
hold on
plot([pi/(T*2*pi) pi/(T*2*pi)],[0 max(abs(FXs2))],'r--','LineWidth',2)
xlabel('f - (Hz)')
set(gca,'FontSize',20)
set(gca,'FontWeight','Bold')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%V=@(z,mu,s) exp(-(z-mu).^2./s);
%s=2000;
V=@(z,mu,s) abs(z-mu)<s;
s=50;
z=1:L;

STEP_fft=2;
Lnow=fix(L/STEP_fft);
%Vn=
%figure
for i=1:Lnow
    Vnow=V(z,STEP_fft*i,s);
    %Vn=Vnow./sum(Vnow);
    Vn=Vnow;
    ys=xs2.*Vn;
    %plot(ys)
    %pause
     FXnow(i,:)=fft(ys,N);
end

figure
ts2=ts(1:STEP_fft:end-1);
pcolor(ts2,f,abs(FXnow)')
shading interp

%%%%%%%%%

xs3=sin(10*ts.^2);
%xs3=sin(0.5*ts.^3);
%ts=T:T:10+T;
%xs3=sin(130./(ts+0.7));

figure 
plot(ts,xs3,'r')

STEP_fft=2;
Lnow=fix(L/STEP_fft);
%Vn=
%figure
for i=1:Lnow
    Vnow=V(z,STEP_fft*i,s);
    %Vn=Vnow./sum(Vnow);
    Vn=Vnow;
    ys=xs3.*Vn;
    %plot(ys)
    %pause
     FXnow(i,:)=fft(ys,N);
end
figure
ts2=ts(1:STEP_fft:end-1);
pcolor(ts2,f,abs(FXnow)')
shading interp