// Blur time envelopes of band-decomposed speech signals
//
// The signal is filtered into several bands and the 
// energy envelope of each one is obtained. They are 
// low-pass filtered to reduce (blur) time details
// and used to reconstruct the signal with blured
// envelopes modulatin noise bands
//
// Requires spectrogram, bands, bandnoise, fftfilter
// 
// --------------------------------
// Federico Miyara
// 2025-09-18


// Clean start
clear
close(winsid())


// Load spectrogram function
exec spectrogram.sci;
// Load band decompose and energy envelope function
exec bands.sci;
// Load band noise generating function
exec bandnoise.sci;
// Load energy envelope function
exec envelope.sci;
// Load white noise generation function
exec whitenoise.sci;
// Load pink noise generation function
exec pinknoise.sci;
// Load brownian noise generation function
exec browniannoise.sci;
// Load mean FFT spectrum function
exec fftmean.sci;
// Load delayed convolution function
exec conv_delay.sci;

// Open wav file containing some speech utterances
path = "R:\Mis documentos\Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\";
caso = 4
select caso
case 1
   fil = path + "tres_frases_consecutivas.wav";
case 2
   fil = path + "Frase (google).wav";
case 3
   fil = path + "Texto (luvvoice).wav";
case 4
   fil = path + "Texto (FSM).wav";
end
[x,Fs] = wavread(fil);
x = x(1,:);
// Time vector
t = (0:length(x)-1)/Fs;
// Signal duration
T = t($)
// RMS value
Xef = sqrt(mean(x.^2))


// FFT window size
N = 4096;
// Mean FFT spectrum of signal
tic
Xfftmean = fftmean(x, Fs, N, "hann", N/2);
toc
// Frequency vector
Fo = (0:N/2-1)*Fs/N;

// Time constant for general envelope extraction
tau = 0.01;
// Compute general envelope with no filtering
tic
Xenv = envelope(x, tau, Fs, 1);
toc
Xenv_ef = sqrt(mean(Xenv.^2))


// Play utterance
playsnd(x, Fs);

// Plot signal and its mean spectrum
scf(1);
clf(1);
gcf().figure_position = [120,130];
subplot(2,1,1)
plot(t,x)
title("Signal")
xlabel("t  (s)"), ylabel("x")
gca().data_bounds = [0 -1; T 1];
subplot(2,1,2)
plot("ll", Fo, Xfftmean)
title("Mean FFT spectrum of signal ")
xlabel("f  (Hz)"), ylabel("Xfftmean")
gca().data_bounds = [20 1e-5;20000 1e-1];

scf(2);
clf(2);
gcf().figure_position = [140,110];
plot(t, x, t, Xenv, "r")
gca().data_bounds = [0 -max(abs(x)); T max(abs(x))];

// Reconstruct signal with envelope and white noise
// Generate white noise
wn = whitenoise(T, Fs);
// Generate white-noise-based signal 
x_env_wn = Xenv.*wn;
x_env_wn = x_env_wn * Xef / stdev(x_env_wn); 
M = max(abs(x_env_wn))
if M>1
   // Prevent signal from being > 1
   x_env_wn = x_env_wn / M;
end
Xef_env_wn = sqrt(mean(x_env_wn.^2));
mprintf("Xef_env_wn = %f\n", Xef_env_wn)
// Play sound 
playsnd(x_env_wn, Fs, 1)

// Reconstruct signal with envelope and pink noise
// Generate pink noise
pn = pinknoise(T, Fs);
// Generate pink-noise-based signal 
x_env_pn = Xenv.*pn;
x_env_pn = x_env_pn * Xef / stdev(x_env_pn);
M = max(abs(x_env_pn))
if M>1
   // Prevent signal from being > 1
   x_env_pn = x_env_pn / M;
end
Xef_env_pn = sqrt(mean(x_env_pn.^2));
mprintf("Xef_env_pn = %f\n", Xef_env_pn)
// Play sound
playsnd(x_env_pn, Fs, 1)

// Reconstruct signal with envelope and brownian noise
// Generate brownian noise
bn = browniannoise(T, Fs);
// Generate brownian-noise-based signal 
x_env_bn = Xenv.*bn;
x_env_bn = x_env_bn * Xef / stdev(x_env_bn);
M = max(abs(x_env_bn))
if M>1
   x_env_bn = x_env_bn / M;
end
Xef_env_bn = sqrt(mean(x_env_bn.^2));
mprintf("Xef_env_bn = %f\n", Xef_env_bn)
// Play sound
playsnd(x_env_bn, Fs, 1)


// BAND TYPES
caso = 5
select caso
case 1
   // Critical bands
   typ = 0;
case 2
   // Equal speech intelligibility 
   typ = -1;
case 3
   // octave 
   typ = 1
case 4
   // One-third octave 
   typ = 3
case 5
   // Two bands
   typ = 20*10^[0 1.5 3]
case 6
   // three bands
   typ = 20*10^[0 1 2 3]
end


// Time constant for computing envelopes
tau = 0.010;
// No low-frequency downsampling of the envelopes
nd = 1;
// Compute band-filtered signals and their envelopes
// fo, fi and fs are column vectors containing the 
// center, lower and upper frequencies of the bands
tic
[y, yenv, fo, fi, fs] = bands(x, typ, Fs, tau, nd);
toc


// RMS value of noise bands
Nef = 1;
// FFT filter window size
N = 4096;
// Random noise
method = "rand";
// Allocate memory
nb = zeros(y);

// Create noise bands
tic
for i = 1:size(nb,1)
   nb(i,:) = bandnoise(fi(i), fs(i), T, Fs, Nef, N, method);
end
toc

// Reconstruct signal with suppressed tonality
// (to certain extent...)
xx = sum(nb.*yenv, 1);
max(xx)
playsnd(xx,Fs,1)

// Signal spectrogram
scf(3);
clf(3);
gcf().figure_position = [150,110]
gcf().figure_size = [1000,600]
subplot(1,2,1)
[z, tto, ffo] = spectrogram(x, 1024, 512, Fs, "hann", "log");
title("Original signal")
gca().data_bounds = [0 0; 5 22000];
xlabel("t  (s)")
ylabel("f  (Hz)")
subplot(1,2,2)
[zz, ttto, fffo] = spectrogram(xx, 1024, 512, Fs, "hann", "log");
title("Signal with spectral blur")
gca().data_bounds = [0 0; 5 22000];
xlabel("t  (s)")
ylabel("f  (Hz)")



// TIME BLUR

// Method 1: Anticipated finite impulse response 
// applied to the envelope

// Delay
d = -5
// Duration of impulse response
L = 2*abs(d) + 1 
// Impulse response
h = [0:abs(d) abs(d)-1:-1:0]/abs(d);

// Initialize smoothed envelopes
yyenv = zeros(size(y,1), size(y,2) + L -1);
// convolve envelopes with bluring impulse response h
tic
for i=1:size(y, 1)
   //a =  conv_delay(yenv(i,:), h, d); 
   yyenv(i,:) = conv_delay(yenv(i,:), h, d); 
end
toc


xxx = sum(nb .* yyenv(:, 1:size(yenv,2)), 1);
max(xxx)
xxx = xxx*Xef/max(xxx);
playsnd(xxx,Fs,1)


// ---------------------------------------
// AGREGAR CONVOLUCIÓN DE LAS ENVOLVENTES
// CON UNA RESPUESTA AL IMPULSO TRIANGULAR
// ANTICIPADA
// ---------------------------------------



// Method 2: 
// Envelope low-pass filtering (time blur)
// Time constant
tau_env = 0.5;
// Laplace transfer function
Hs = 1/(1 + tau_env*%s);
// Bilinear transformation
s = 2*Fs*(%z - 1)/(%z + 1);
// Z transfer function
Hz = horner(Hs,s)
NUM = Hz.num
DEN = Hz.den
// Initialize smoothed envelopes
yyenv = zeros(y);
// Filter envelopes
tic
for i=1:size(y, 1)
   yyenv(i,:) = filter(NUM, DEN, yenv(i,:)); 
end
toc

// Reconstruct signal with suppressed tonality
x_blur = sum(nb.*yyenv, 1);
max(x_blur)
playsnd(x_blur,Fs,1)


// Plot critical-band envelopes
M = max(abs(y)); 
scf(4);
clf(4);
m = 6;
n = 4
for i = 1:m*n
   // Modified index to plot subplots by columns
   // instead of rows
   I = 1 + floor((i-1)/n) + m*modulo(i-1, n);
   subplot(6,4,i)
   plot(t, yenv(I,:), "r")
   //plot(t, y(I,:))
   ylabel("y")
   legend(string(floor(fo(I))) + " Hz")
   gca().data_bounds = [0, -M ; max(t) M];
end

// Plot original signal and signal reconstructed
// suppressing tonality
scf(2);
clf(2);
subplot(2,1,1)
plot(t,x)
subplot(2,1,2)
plot(t,xxx)

