% espectrograma 

close all
clear all

% parametros
N = 512;
fs = 8000;
M = 20;
T = N*M;
t = ((1:T)-1)/fs;

win = hanning(T);

% seniales
f1 = 750;  % Hz
f2 = 2500; % Hz

ruido = 2*rand(T,1)-1;
senial1 = sin(2*pi*f1*t(:,1:T/2)');
senial2 = sin(2*pi*f2*t(:,1:T/2)');
senial = [senial1; senial2] / 2; 
senial = senial + 0.01*ruido;

wavwrite(senial,fs,'senial2.wav');

% espectro de la senial 
X = 4/T * fft(senial.*win);

figure;
subplot(211)
plot(t,senial,'k')
title('Senial')
xlabel('Tiempo (s)')
ylabel('Amplitud')
axis([0 t(end) -1 1])
subplot(212)
plot(0:fs/T:fs/2,20*log10(abs(X(1:T/2+1))),'k')
xlabel('Frecuencia (Hz)')
ylabel('Magnitud (dB)')

print -depsc 'senial_no_estacionaria.eps'

disp('pulse una tecla para seguir ...')
pause

[S, f, t] = specgram(senial, N, fs, N, 0);

figure;
imagesc(t,f,20*log10(abs(S)))
colormap(abs(1-gray))
axis xy
title('Espectrograma')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

print -depsc 'espectrograma.eps'


%% espectrograma de banda angosta y de banda ancha

[y, fs] = wavread('f24cb9_4.wav');

N1 = 2048;
N2 =  256;

[S1, f1, t1] = specgram(y, N1, fs, hanning(N1), N1-N1/4);
[S2, f2, t2] = specgram(y, N2, fs, hanning(N2), N2-N2/4);

fmax = 10000;
ind1_max = find(f1 > fmax, 1, 'first') - 1;
ind2_max = find(f2 > fmax, 1, 'first') - 1;

figure;
subplot(211)
imagesc(t1,f1(1:ind1_max),20*log10(abs(S1(1:ind1_max,:))))
colormap(abs(1-gray))
axis xy
title('Espectrograma de banda angosta')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

%figure;
subplot(212)
imagesc(t2,f2(1:ind2_max),20*log10(abs(S2(1:ind2_max,:))))
%imagesc(t2,f2,20*log10(abs(S2)))
colormap(abs(1-gray))
axis xy
title('Espectrograma de banda ancha')
xlabel('Tiempo (s)')
ylabel('Frecuencia (Hz)')

print -depsc 'espectrograma_banda_angosta_banda_ancha.eps'
