// Testing function spectral_blur
//
// --------------------------------
// Federico Miyara
// 2025-05-29

// Clean start
clear
close(winsid())

caso = 2
select caso
case 1
   // Open wav file containing some speech utterances
   path = "R:\Mis documentos\Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\"
   name = "tres_frases_consecutivas.wav"
   fil = path + name
   nameblur = part(name,1:length(name)-4) + "_blur.wav"
   [x, Fs] = wavread(fil);
   // Keep only left channel
   x = x(1,:);
   t = [0:length(x)-1]/Fs;
case 2
   path = "R:\Mis documentos\Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\"
   name = "Enigma (FSM).wav"
   fil = path + name
   nameblur = part(name,1:length(name)-4) + "_blur.wav"
   [x, Fs] = wavread(fil);
   // Keep only left channel and from t = 1 s to t = 4.5 s
   x = x(1, Fs:Fs*4.5);
   t = [0:length(x)-1]/Fs;   
case 3
   Fs = 44100;
   fo = 1000
   T = 3
   t = 0:1/Fs:T;
   x = zeros(t);
   x(Fs+1:2*Fs) = 0.3*sin(2*%pi*fo*(2+t(1:Fs)));  
end

// Get original signal spectrogram
exec spectrogram.sci;
N = 1024;
M = N/2;
w = "hann";
grph = "none"; 
cmap = "jet";
DR = 100;
[z, t1, f] = spectrogram(x, N, M, Fs, w, grph, cmap, DR);

// Length of FFT analysis window
N = 1024;
// Critical band widening type                 
typ = "crit"; 
// Width of spectral line widening window expressed
// as a fraction of a critical band          
L = 1;                   
// Load spectral_blur function
exec spectral_blur.sci;
// Get spectrally blurred signal and its spectrogram
tic
[xs, zs, z, f, ts] = spectral_blur(x, Fs, N, typ, L);
toc

// Play original and spectrally blurred signal
playsnd(x,Fs,1)
playsnd(xs,Fs,1)

// Save blurred signal
wavwrite(xs, Fs, path + nameblur)
 
scf(1);
clf(1);
gcf().figure_position = [100,75];
gcf().color_map = jet(64);
gcf().background = -2;
Sgrayplot(t1, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
title("Spectrogram of original signal")
xlabel("t  (s)")
ylabel("f  (Hz)")

scf(2);
clf(2);
gcf().figure_position = [711,75];
gcf().color_map = jet(64);
gcf().background = -2;
Sgrayplot(t1, f, 20*log10(zs'), zminmax=[-DR 0], colout=[-1 -1]);
title("Spectrogram of spectrally-smoothed signal")
xlabel("t  (s)")
ylabel("f  (Hz)")

// Spectrograms to print
scf(3);
clf(3);
gcf().figure_position = [500,75];
gcf().figure_size = [650,400];
gcf().color_map = jet(64);
gcf().background = -2;
subplot(1,2,1)
  Sgrayplot(t1, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
  // title("Spectrogram of original signal")
  xlabel("t  (s)")
  ylabel("f  (Hz)")
subplot(1,2,2)
  Sgrayplot(t1, f, 20*log10(zs'), zminmax=[-DR 0], colout=[-1 -1]);
  // title("Spectrogram of spectrally-smoothed signal")
  xlabel("t  (s)")
  ylabel("f  (Hz)")

MAX = max(max(abs(x)), max(abs(xs)))
scf(4);
clf(4);
gcf().figure_position = [500,55];
subplot(2,1,1)
plot(t, x)
xlabel("t  (s)")
ylabel("y_original")
gca().data_bounds = [0 -MAX; ts(length(t)) MAX];
subplot(2,1,2)
plot(ts(1:length(t)), xs(1:length(t)))
xlabel("t  (s)")
ylabel("y_smoothed")
gca().data_bounds = [0 -MAX; ts(length(t)) MAX];
