function [xs, zs, z, f, ts] = spectral_blur(x, Fs, N, L)
   // Smooth the spectrum of a signal
   //
   // This function gets a signal with similar time
   // envelope to a given signal but widened spectral 
   // lines, which results in a smoothed or blurred
   // spectrogram. 
   // This may be used to keep the general spectral 
   // distribution while reducing the intelligibility
   // or spectral identity of a signal, for instance
   // a speech signal 
   // 
   // Usage:
   //       [xs, zs, z, f, ts] = spectral_blur(x, Fs, N, L)
   // where
   //       x:  Signal vector
   //       Fs: Sample rate (Hz)
   //       N:  Length of analysis window
   //       L:  Length of spectral line widening
   //       xs: Spectrally-smoothed signal vector
   //       zs: Spectrally-smoothed spectrogram
   //       z:  Original spectrogram
   //       f:  Frequency bin vector (Hz)     
   //       ts:  Time slices vector (s)
   //
   // First the real spectrogram is obtained, which
   // consists of a series of spectra, each one 
   // corresponding to a time slice. Second, the 
   // spectrum at each time slice is convolved with
   // a Hann window of length L. The result is a
   // smoothed spectrum. Third, the inverse FFT of
   // each smoothed spectrum (with a random phase 
   // added) is computed, yielding a series of 
   // short overlapping signals. An overlapp-add 
   // procedure is applied to get the spectrally 
   // smoothed signal window.
   //
   // -------------------------------------
   // Author: Federico Miyara
   // Date:   2025-03-20  
   
   rhs = argn(2)
   if rhs<4
      L = 16
   end
   if rhs<3
      N = 1024
   end
   if rhs<2
      Fs = 44100
   end
   
   // For testing purposes only --REMOVE--
   if 1==2
      
   caso = 2
   select caso
   case 1
      // Open wav file containing some speech utterances
      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
   
   end

   // Time vector
   tt = (1:length(x))/Fs;
   
   // ----------------------------------------- 
   // 1 ORIGINAL SPECTROGRAM
   // -----------------------------------------

   // Load spectrogram function
   exec spectrogram.sci;

   // Spectrogram parameters
   // Overlapping
   M = N/2
   // FFT window
   w = "hann" 
   // Type of amplitude representation
   grph = "none"

   // Compute spectrogram 
   [z, t1, f] = spectrogram(x, N, M, Fs, w, grph);

   // Compute the absolute maximum and minimum
   [zmax, I] = max(z)
   [zmin, J] = min(z)

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

   // 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 
   for k=1:size(z,2)
      yy(:,k) = conv(z(:,k), W);
   end
   // Compensate offset due to the widening
   // window and truncate length to N/2 
   y = yy(L/2+1:L/2+N/2,:);

   // 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;


   // ----------------------------------------- 
   // 3 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;

endfunction
