// Several band noises of different bandwidths and
// the same level. The noises are plotted in the
// time domain, in the frequency (spectral) domain,
// and listened to,
//
// -----------------------------------
// Federico Miyara
// 2024-11-03

clear, clc, close(winsid())

// Load functions
exec bandnoise.sci;
exec zbark.sci;
exec zkrab.sci;

// Band center frequency
fo = 800
// Band center critical band rate
zo = zbark(fo)
// Frequency limits of critical band centered at fo
fo1 = zkrab(zo - 0.5)
fo2 = zkrab(zo + 0.5) 
// Duration
T = 3
// Sample rate
Fs = 44100;
// Effective value
Xef = 0.2;
// Length of filtering window
N = 8192;

// Crítical band proportion vector
fracBC = [0.2 0.4 0.7 1.0 1.5 2];

// Bandwidth vector in Hz
df = (fo2 - fo1)*fracBC;

// Create band noises
// If method = "perf" then spectrally perfect 
// (flat) noise is created. Otherwise it will
// be flat only in a long-term average fashion
// Note that spectrally perfect noise will
// present audible repetition pattern unless
// n >= 65536 
method = "perf";
for i=1:length(df)
   f1 = fo - df(i)/2;
   f2 = fo + df(i)/2;
   x(i,:) = bandnoise(f1, f2, T, Fs, Xef, N, method);
end

// Frequency vector
f = (0:N/2-1)*Fs/N;
// Spectra of signals
X = abs(fft(x(:,1:N), -1, 2))(:,1:N/2)*2/N;

// Create plots of the first 100 ms of the six
// signals
scf(1);
clf(1);
TT = 0.1
t = 0:1/Fs:TT;
for i=1:length(df)
   subplot(length(df)/2,2,i)
   if fracBC(i)==1
      colr = "r";
   else
      colr = "b";
   end 
   plot(t, x(i,1:length(t)), colr);
   legend(string(fracBC(i))+" BC", 2)
   gca().data_bounds = [0, -max(x); TT, max(x)]
end
 
// Create plots of the spectra of the six signals
// around the center frequency
clf(2);
for i=1:length(df)
   subplot(length(df)/2,2,i)
   if fracBC(i)==1
      colr = "r";
   else
      colr = "b";
   end 
   plot("ln", f, X(i,:),colr);
   legend(string(fracBC(i))+" BC", 2)
   gca().data_bounds = [fo/1.5, 0; fo*1.5, max(X)];
end

// Check effective values
Xef = sqrt(mean(x.^2,"c"));
mprintf("Xef(%i) = %.2f \n", (1:6)', Xef)

// Play noises from narrow to wide. The fourth
// noise is a critical-band-wide one
for i=1:length(df)
   playsnd(x(i,:), Fs, 1)
end

// Save noises as WAV
// (To really save replace path and change "n" to "y")
S = "n";
if S=="y"
   path = "R:\Mis documentos\Acústica\Cursos\Psicoacústica\Sonidos\";
   name = "ruido_fo_800Hz_" + string(fracBC) + "BC.wav"; 
   fil = path + name;
   for i=1:length(df)
      wavwrite(x(i,:), Fs, fil(i));
   end
end
