
// Get a signal with similar time envelope to a given 
// signal but widened spectral lines, which results 
// in a smoothed or blurred spectrogram.
//
// First the real spectrogram is obtained, resulting
// in a series of spectra, each one corresponding to 
// a time slice. Second, the spectrum at each time 
// slice is convolved with a Hann window, resulting
// in a smoothed spectrum. Third, the inverse FFT of
// each smoothed spectrum to which a random phase is
// attached is computed, yielding a series of short
// overlapping signals. An overlapp-add procedure is
// applied to get the spectrally smoothed signal    
// window at each time where the 
//
// x:  Original signal
// z:  Spectrogram of the original signal
// W:  Smoothing window
// y:  Smoothed spectrogram
// yo: Normalized smoothed spectrum to get equal
//     RMS as the original spectrum 
// xb: Spectrally blurred, noise-like signal
//
// --------------------------------
// Federico Miyara
// 2025-03-16


// Clean start
clear
close(winsid())

// Load spectrogram function
exec spectrogram.sci;

// Open wav file containing some speech utterances
caso = 2
select caso
   case 1
      path = "R:\Mis documentos\Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\"
      fil = path + "tres_frases_consecutivas.wav"
      [x,Fs] = wavread(fil);
      x = x(1,:);
   case 2
      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

// Time vector
t = (1:length(x))/Fs;

// Optionally, play the sound
play = "n"
if play=="y"
   playsnd(x,Fs)
end


// ----------------------------------------- 
// 1 ORIGINAL SPECTROGRAM
// -----------------------------------------

// Spectrogram parameters
// Size of FFT window
N = 1024
// Overlapping
M = N/2
// FFT window
w = "hann" 
// Type of amplitude representation
grph = "log"
// Color map
cmap = "jet"

// Create figure 
scf(1);
clf(1);
gcf().figure_position = [100,75];

// Draw spectrogram and get spectral data
[z, t1, f] = spectrogram(x, N, M, Fs, w, grph, cmap);
title("Spectrogram of original signal")
xlabel("t  (s)")
ylabel("f  (Hz)")

// Compute and place on the spectrogram the 
// absolute maximum and minimum
[zmax, I] = max(z)
[zmin, J] = min(z)
plot(t1(I(2)), f(I(1)), "o")
plot(t1(J(2)), f(J(1)), "o")


// ----------------------------------------- 
// 2 SMOOTHED SPECTROGRAM
// -----------------------------------------

// Spectral line widening size
L = 32
// Symmetric Hann window for widening
// Might be other type, such as Gaussian,
// and its size could be also made a
// function of frequency, por instance 
// to match the ear frequency resolution
W = window("hn", L);

// Widening window mean and RMS 
Wmean = mean(W) 
WRMS = sqrt(mean(W.^2))

// Widening window ideal mean and RMS values
Wmean_ideal = 0.5
WRMS_ideal =sqrt(0.5^2 + 0.5^2/2)

// Widen lines convolving with symmetrical 
// widening window 
clear yy
tic
for k=1:size(z,2)
   yy(:,k) = conv(z(:,k), W);
end
toc
// Compensate offset due to the widening
// window and truncate length to N/2 
y = yy(L/2+1:L/2+N/2,:);

// Limit dynamic range to remove graphic 
// issues and problems with log(0)
DR = 100
y = max(y, max(y)*10^(-DR/20));

// RMS of original and blurred spectrogram
zRMS = sqrt(mean(z.^2))
yRMS = sqrt(mean(y.^2))

// Ratio between blurred and original RMS 
KRMS = yRMS/zRMS

// Normalize the smoothed spectrogram to the 
// RMS value of the original spectrogram
yo = y/KRMS;

// Time vector
t2 = (0:size(y,2)-1)*M/Fs;

// In order to plot yo with the same color
// reference as z, the value of yo at the 
// (t, f) point in the spectrogram were z 
// was maximum will be replaced by zmax,
// creating a new y1 spectrogram. Using
// the same range of values will ensure that
// the color reference will be the same.
// A single more intense point is barely
// noticeable 
y1 = yo;
//y1(I(1),I(2)) = zmax;
//y1 = max(y1, max(y1)*10^(-DR/20)); 

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


// ----------------------------------------- 
// 3 COMPARISON BETWEEN ORIGINAL AND
//   SMOOTHED SPECTRA
// -----------------------------------------

// Index of a specific time (where the spectrum
// reaches its maximum) 
m = I(2)

// Compare original and blurred spectra at 
// the specific time of the maximum spectrum  
scf(3);
clf(3);
gcf().figure_position = [400,55];
subplot(2,1,1)
plot((0:N/2-1)*Fs/N, z(:,m))
gca().data_bounds = [0 0; 23000 0.4];
xlabel("f  (Hz)")
ylabel("Original spectrum")
subplot(2,1,2)
plot((0:N/2-1)*Fs/N, yo(:,m))
gca().data_bounds = [0 0; 23000 0.4];
xlabel("f  (Hz)")
ylabel("Smoothed spectrum")

// Verify RMS value of one particular
// frame where the spectrum reaches its 
// absolute maximum
zoo = N/2 * ifft([z(:, m); 0; z($:-1:2, m)]);
sqrt(mean(zoo.^2))
yoo = N/2 * ifft([yo(:, m); 0; yo($:-1:2, m)]); 
sqrt(mean(yoo.^2))



// ----------------------------------------- 
// 4 CONVERT UNILATERAL SMOOTHED SPECTRA
//   TO FFT FOR EACH TIME SLICE
// -----------------------------------------

// Assign a random phase and convert the 
// unilateral spectrum that corresponds to 
// each time slice into a complex bilateral,
// conjugate-symmetrical spectrum
Xo = zeros(N, size(yo, 2));
for k = 1:size(yo,2)
   Xo(1:N/2, k) = yo(:, k) .* exp(%i*2*%pi*rand(N/2, 1));
   Xo(1, k) = 0;
   Xo(N/2+1, k) = 0;
   Xo($:-1:$-N/2+2, k) = conj(Xo(2:N/2, k));
end

// Get signals corresponding to each time 
// slice. Since the spectrogram removes
// the N/2 factor implicit in the FFT to
// get true spectral compenent amplitude,
// this factor must be applied to the 
// inverse FFT  
xo = N/2 * fft(Xo, 1, 1);

// Asymmetric Hann window to cross-fade
// between consecutive signals
h = window("hn", N+1)(1:N)';
hh = repmat(h,1,size(xo,2));
// Replace first and last halves bt ones
hh(1:N/2,1) = ones(N/2,1);
hh(N/2+1:$,size(xo,2)) = ones(N/2,1);

// Multiply each signal by the respective
// window
xo = xo.*hh;

// Concatenate odd and even signals  
xodd = [];
xeven = zeros(N/2,1);
for k=1:2:size(xo,2)
   xodd = [xodd; xo(:,k)];
   xeven = [xeven; xo(:,k)];
end

if length(xeven)<length(xodd)
   xeven = [xeven; zeros(N/2,1)];
elseif length(xodd)<length(xeven)
   xodd = [xodd; zeros(N/2,1)];
end

xs = xeven + xodd;
ts = (0:length(xs)-1)/Fs;

scf(4);
clf(4);
gcf().figure_position = [500,55];
subplot(2,1,1)
plot(t, x)
xlabel("t  (s)")
ylabel("y_original")
subplot(2,1,2)
plot(ts, xs)
xlabel("t  (s)")
ylabel("y_smoothed")

playsnd(x, Fs, 1)
playsnd(xs, Fs)
playsnd(max(abs(x))/max(abs(xs))* xs, Fs)

      
   


//////////////////////////////////////////////
//
// VER por qué la amplitud se reduce tanto
// VER por qué hay un retardo de 512 muestras
//
//////////////////////////////////////////////


