
// Prueba de la función notched noise

// Tasa de muestreo
Fs = 44100
// Duración
T = 20
// Frecuencia central
f = 1000
// Ancho de banda
BW = 100
// Ventana de análisis
N = 32768

// Cantidad de ventanas completas
m = floor(T*Fs/N)

// Cargar la función
exec notched_noise.sci;

// Generación de ruido blanco espectralmente
// perfecto
yp = notched_noise(f, BW, T, Fs, "p");

// Generación de ruido blanco aleatorio
y = notched_noise(f, BW, T, Fs);


// Reproducción del ruido
playsnd(0.2*yp(1:100000), Fs, 1);
playsnd(0.2*y(1:100000), Fs, 1); 

// Cálculo del espectro promedio de todas 
// las ventanas completas (ambos ruidos)
exec fftmean.sci;
[Yp2, Yp, ypmin, ypmax, fo] = fftmean(yp(1:m*N), Fs, N, "box", N/2);
[Y2, Y, ymin, ymax, fo] = fftmean(y(1:m*N), Fs, N, "box", N/2);

// Cálculo del espectro de una sola ventana
// (ambos ruidos) 
Ypo = abs(fft(yp(1:N))(1:N/2))/sqrt(2);
Ypo(1) = Ypo(1)/N;
Ypo(2:$) = 2*Ypo(2:$)/N;

Yo = abs(fft(y(1:N))(1:N/2))/sqrt(2);
Yo(1) = Yo(1)/N;
Yo(2:$) = 2*Yo(2:$)/N;


scf(1);
clf(1);
gcf().figure_size = [626,416];  
gcf().figure_position = [200,100];
plot("nn", fo, 20*log10(Yp2), fo, 20*log10(Ypo)); 
title("Ruido espectralmente perfecto")
xlabel("f   (Hz)")
ylabel("Yp   (dB)")
gca().data_bounds = [   800, max(min(20*log10(Yp2)) - 5, -360) 
                     1200, max(20*log10(Yp2)) + 5];

scf(2);
clf(2);
gcf().figure_size = [626,416];
gcf().figure_position = [220,80]
plot("nn", fo, 20*log10(Y2), fo, 20*log10(Yo)); 
title("Ruido blanco aleatorio")
xlabel("f   (Hz)")
ylabel("Y   (dB)")
gca().data_bounds = [   800, max(min(20*log10(Y2)) - 5, -360) 
                     1200, max(20*log10(Y2)) + 5];
