
// Prueba de la función trackfilter
//
// Esta función filtra una señal siguiendo la evolución
// de una frecuencia dada como dato en función del tiempo.
// El filtro puede ser pasabajos, pasaaltos, pasabanda, 
// rechazabanda, y pasabanda o rechazabanda de múltiples
// armónicos
//
// ----------------------------
// Federico Miyara
// 2008-03-04
// 2021-08-29
// 2023-12-04
// 2025-09-27

clear
close(winsid())

// Tamaño de la FFT
N = 2048;

// Tasa de muestreo
Fs = 44100;

// Duración
T1 = 4;

// Vector de tiempos
tt = [0:1/Fs:T1];

// Señal de ruido blanco
xx = rand(tt, "normal");
xx = 0.5*xx/max(abs(xx));

// Gráfica de la señal
scf(1);
clf(1);
gcf().figure_position = [200 120];
plot(tt, xx);
xlabel('t  [s]')
ylabel('x')

// Audición de la señal
// playsnd(x, Fs, 1)

// Tipo de acción filtrante
mode = "bp5";
//mode = "lp";

// Ancho de banda dado como porcentaje de la frecuencia
// del filtro si aparece '%' en mode o como porcentaje
// de la tasa de muestreo si no aparece 
B = 1;

// Evolución de la frecuencia fundamental
// Cantidad de muestras de frecuencia 
L = 100;

// Vector de tiempos discretos
T = [0:L-1]*T1*Fs/(L-1);

// Frecuencia máxima
Fmax = 5000;

// Variable proporcional a T que fluctúa de 0 a 1
u = [0:L-1]/(L-1);

// Vector de frecuencias
F = (u + sin(u*2*2*%pi)/2/%pi) * Fmax/Fs;

// Gráfica de la frecuencia fundamental
scf(2)
clf(2);
gcf().figure_position = [220 100];
plot(T/Fs, F*Fs)
xlabel('t  [s]')
ylabel('f  [Hz]')

// Filtrado de la señal
exec trackfilter.sci;
x = real(trackfilter(xx, F, T, N, mode, B));
Xmax = max(abs(x));
if Xmax > 1
   x = 0.5 * x/Xmax;
end

// Gráfica de la señal filtrada
scf(3);
clf(3);
gcf().figure_position = [240 80];
plot([0:length(x)-1]/Fs, x)
title('Señal filtrada')
xlabel('t  [s]')
ylabel('y')

// Espectrograma de la señal filtrada
exec spectrogram.sci;
N = 1024;
Fs = 44100;
M = 512;
w = "hann";
grph = "logcmap";
[z, t, f] = spectrogram(x, N, M, Fs, "hann", "logcmap");

// Graficas de la señal, la frecuencia a seguir
// y el espectrograma de la señal filtrada
scf(4);
clf(4);
gcf().figure_position = [260 60];
subplot(3,1,1)
plot(tt, xx);
xlabel('t  [s]')
ylabel('x')
subplot(3,1,2)
plot(T/Fs, F*Fs)
xlabel('t  [s]')
ylabel('f  [Hz]')
subplot(3,1,3)
cmap = jet(64)
gcf().color_map = cmap;
gcf().background = -2;
Sgrayplot(t, f, 20*log10(z'));
 
            
// Audición de la señal filtrada
playsnd(0.2*x/max(abs(x)),Fs,1)

// Guardado de la señal filtrada
//wavwrite(x, Fs, 'e:\ondas\pruebatrackfilter.wav')

