clear;
% Dos señales sinusoidales de ~50 Hz

% Queremos una resolución frecuencial de 0.1 Hz
df=0.1;
D=1/df;  

% La frecuencia de muestreo tiene que ser > 2*50
fm=200;ts=1/fm;  
N=D/ts;  				
N=2^(ceil(log(N)/log(2)));     % aproximo a la potencia de 2 más cercana
ts=D/N;   				% recalculo ts y fm
fm=1/ts;
fprintf(1,'\n***********************************************************************\n');
fprintf(1,'Caso 1: Hago la FFT de todos los puntos de señal en %.2f segundos\n',D);
fprintf(1,'La frecuencia de muestreo es de %.2f Hz\n',fm);
fprintf(1,'El número de puntos de la FFT es %d\n\n\n',N);

t=0:ts:D-ts/2;

x=sin(2*pi*49.5*t);
% 0.3*sin(2*pi*150*t)+0.1*sin(2*pi*300*t);

X=fft(x);
X=fftshift(X);
f=-N/2:0;
f=[f 1:N/2-1];f=f*df;

stem(f,abs(X));

disp('Presione una tecla para continuar');
pause;

% Llenado de ceros
% La resolución es 0.1 Hz
% Muestreo durante 10 ciclos de la señal de 50Hz = 200ms
D=200e-3;

% Tomamos como criterio que 8M es ~Nz
% M es ~50/df --> 8*50/df = Nz
% Para que Nz sea potencia de 2 --> n=ceil(log(Nz)/log(2)) --> Nz=2^n;


Nz = 8*50/df;
n = ceil(log(Nz)/log(2));Nz=2^n;

ts=1/(df*Nz);fm=1/ts;   % recalculo los valores de ts y fm

fprintf(1,'\n***********************************************************************\n');
fprintf(1,'Caso 1: Determino el espectro de la señal mediante el llenado de ceros\n');
fprintf(1,'Muestreo la señal durante %.2f segundos (~%.2f ciclos)\n',D,D*50);
fprintf(1,'La frecuencia de muestreo es de %.2f Hz\n',fm);
fprintf(1,'El número de puntos de la FFT es %d\n\n\n',Nz);

t=0:ts:D-ts/2;
x=sin(2*pi*49.5*t);
w=hanning(length(t))';

x=x.*w;
X=fft(x,Nz);
X=fftshift(X);
Xmax=max(abs(X));
f=[-Nz/2:0 1:Nz/2-1];f=f*df;

plot(f,20*log10(abs(X/Xmax)));grid;zoom;