
// Prueba de la función morphfilter
//
// Esta función filtra una señal con una respuesta en frecuencia
// variable en el tiempo según una definición dada por una matriz
// cuyas columnas representan las respuestas en frecuencia 
// correspondientes a sucesivos instantes de tiempo. Los instantes
// en los que se especiufica cada respuesta en frecuencia se dan 
// en un vector aparte o de lo contrario se los supone equiespaciados
// a lo largo de la señal.
//
// ----------------------
// Federico Miyara
// 2021-07-21
// 2023-10-30


clear

// Duración de la señal [s]
TT = 3;

// Tasa de muestreo [Hz]
Fs = 44100;

// Tamaño de las FFT
N = 4096;

// Vector de tiempos
t = [0:1/Fs:TT];

// Ruido blanco
x = rand(t,"n");

// Longitud original de la señal 
Lo = length(x);

// Vector de índices (tiempos discretos) de las muestras
// de respuesta en frecuencia 
M = 100;
T = (Lo + 2048)*[1:M]/M;

// Especificación de las ventanas filtrantes. Se dan M 
// ventanas pasabanda de ancho de banda relativo constante,
// equiespaciadas en el tiempo y con frecuencias superior e
// inferior de corte crecientes con el tiempo.
M = 100;
for k=1:M
   Fsup = round(k/M*N/3);
   Finf = round(Fsup*3/4);
   W(:,k) = [zeros(Finf, 1); ones(Fsup-Finf, 1); zeros(N/2 - Fsup, 1)];
end

// Señal filtrada
exec morphfilter.sci;
z = real(morphfilter(x,W,T,N));

// Prueba de sonido de la señal de entrada
playsnd(0.1*x/max(abs(x)), 44100)

// Prueba de sonido de la señal filtrada
playsnd(0.1*z/max(abs(z)),44100)

// Escritura de la señal filtrada como archivo wav 
// wavwrite(0.9*z/max(abs(z)), Fs, 'e:\Ondas\morphfilter2.wav')

// Gráfica en 3D de la respuesta en frecuencia dependiente del tiempo
scf(1);
clf(1);
gcf().figure_position = [200, 110];
mesh(1:N/2, 1:5:M,  W(:,1:5:M)')

// Espectrogramas de la señal de entrada y de la señal filtrada 
exec spectrogram.sci;
scf(2);
clf(2);
gcf().figure_position = [220, 90];
subplot(2,1,1)
spectrogram(x, 1024, 512, Fs, "hann", "log");

subplot(2,1,2)
spectrogram(z, 1024, 512, Fs, "hann", "log");

// Oscilograma de la señal filtrada
scf(3);
clf(3);
gcf().figure_position = [240, 70];
plot([0:length(z)-1]/Fs, real(z)) 
gca().grid = [1 1]

