%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% EFFECT OF WINDOW SHAPE ON THE SPECTRUM OF SINUSOIDAL SIGNALS
%% (Dimitris, page 408)
%%
%% a) La ventana rectangular tiene la MAYOR resolucion pero el MAYOR
%% ripple.
%%
%% El resto de las ventanas (Bartlett,Hann y Hamming) tienen MENOR
%% ripple a expensas de tambien MENOR resolucion.
%%
%% b) La ventana de Bartlett (o triangular) tiene MAYOR ripple respecto
%% de Hann y Hamming pero mejor resolucion que Hamming.
%%
%% c) La ventana de Hann tiene el MENOR ripple pero una resolucion
%% levemente menor que Bartlett y Hamming.
%%
%% d) La ventana de Hamming tiene la MAYOR resolucion, a expensas de un
%% mayor ripple respecto de la ventana de Hann.
%%
%%
%% La ventada de observacion de largo T0 debe satisfacer:
%%
%% T0 > 1/(f3-f2)
%%
%% para la ventana rectangular para poder resolver las frecuencias f2 y
%% f3 y
%%
%% T0 > 2 1/(f3-f2)
%%
%%  para el resto de las ventanas (ya que el lobulo principal tiene un
%%  ancho que es aproximadamente el doble que el de la ventana
%%  rectangular).
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clear all
close all


% SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f1  =  1  % Hz
f2  =  3
f3  =  4

A1=1.0;
A2=1.0;
A3=1.0;

fN=10*f3  % donde quiero que este Nyquist

dt=1.0/(2*fN);  %fS=2*fN



%% Ventana de observacion

%% T0 > 1/(f3-f2) para rectangular y T0 > 2 1/(f3-f2) para el resto

T0=1/(f3-f2);     % Ventana de observacion minima para obtener
		  % resolucion en la ventana rectangular.(Este valor es
		  % chico para las ventanas distintas a la rectangular)

%T0=2/(f3-f2);    % Ventana de observacion minima para obtener resolucion en
		  % todas las ventantas


k=2; % aumenta k, aumenta el largo de la ventana de observacion =>
				% aumenta la resolucion del espectro
T0=k*T0

t=(0:dt:T0); % intervalo de observacion (ventana rectangular)


% SENAL

s1=A1*sin(2*pi*f1*t);
s2=A2*sin(2*pi*f2*t);
s3=A3*sin(2*pi*f3*t);

s=s1+s2+s3;

N=length(s);


% VENTANA SOBRE LA SENAL

Nw=round(T0/dt);  % numero de coeficientes de la ventana

if(mod(Nw,2) == 0)
  Nw=Nw+1;     % numero impar de coeficientes de la ventana
endif
Nw
Mw=(Nw-1)/2    % mitad del numero de coeficientes de la ventana

sbar=s; shann=s; shamm=s;

window=bartlett(Nw);

for i=1:Mw
  sbar(i)=s(i)*window(i);
  sbar(N-Mw-1+i)=s(N-Mw-1+i)*window(Mw-1+i);
     endfor

window=hanning(Nw);
for i=1:Mw
  shann(i)=s(i)*window(i);
       shann(N-Mw-1+i)=s(N-Mw-1+i)*window(Mw-1+i);
     endfor

window=hamming(Nw);
for i=1:Mw
  shamm(i)=s(i)*window(i);
       shamm(N-Mw-1+i)=s(N-Mw-1+i)*window(Mw-1+i);
     endfor


% ESPECTROS DE AMPLITUD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Npad=2**10; % zero padding

spad      = s   ;      spad(Npad) = 0;      Spad = abs(fft(spad));
spadbar   = sbar;   spadbar(Npad) = 0;   SpadBar = abs(fft(spadbar));
spadhann  = shann; spadhann(Npad) = 0;  SpadHann = abs(fft(spadhann));
spadhamm  = shamm; spadhamm(Npad) = 0;  SpadHamm = abs(fft(spadhamm));

fS=1/dt;
f=(0:Npad-1)*(fS  /Npad);

final = 10; % fS/2.0;       % frecuencia final para plotear  (Hz) [0,fS)

norm=2/N; % para ver las amplitudes de los senos de la senal (el dos es
	 % porque junto la energia despues de Nyquist en el intervalo [0,Nyquist]


% GRAFICOS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)

subplot(2,2,1)
plot(t,s,'k');title('w_n, Rectangular')
ylabel('s_n x w_n')
subplot(2,2,2)
plot(t,sbar,'r');title('w_n, Bartlett')
subplot(2,2,3)
plot(t,shann,'g');title('w_n, Hanning')
xlabel('tiempo (s)')
ylabel('s_n x w_n')
subplot(2,2,4)
plot(t,shamm,'b');title('w_n, Hamming')
xlabel('tiempo (s)')

figure(2)

subplot(2,2,1)
plot(f,norm*Spad,'k');title('Rectangular')
ylabel('2/N |Sw|')
axis([0  10 -0.05 1.2])
subplot(2,2,2)
plot(f,norm*SpadBar,'r');title('Bartlett')
axis([0  10 -0.05 0.8])
subplot(2,2,3)
plot(f,norm*SpadHann,'g');title('Hanning')
xlabel('frecuencia (Hz)')
ylabel('2/N |Sw|')
axis([0  10 -0.05 0.8])
subplot(2,2,4)
plot(f,norm*SpadHamm,'b');title('Hamming')
xlabel('frecuencia (Hz)')
axis([0  10 -0.05 0.8])


% OUTPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

