function [Xrho, fo, X, Xenv, tenv] = envelope_corr(x, Fs, typ, T, Fsenv)
   // Correlation between spectral envelopes of N signals
   //
   // Usage:
   //       [Xrho, fo, X, Xenv, tenv] = envelope_corr(x,Fs,typ,T,Fsenv)
   // where
   //       x:     1-D Cell array containing signal vectors 
   //              or paths to wav files (can be mixed). In
   //              this case only the left channel will be
   //              used for stereo files 
   //       Fs:    Vector of the same size as the cell x
   //              so that Fs(n) is the sample rate of x{n}
   //       typ:   Spectral band type:
   //              equal speech intelligibility (-1), 
   //              critical (0), octave (1), one-third 
   //              octave (3), etc.; also any custom set 
   //              of frequencies which define adjacent 
   //              frequency bands, arranged as a vector.
   //              If not in increasing order, the vector
   //              will be reordered. 
   //       T:     Duration in s of the relevant part of
   //              the signals to analyze
   //       Fsenv: Downsampled envelope sample rate
   //       Xrho:  2-D cell array of vectors such that
   //              Xrho{n,m}(k) is the maximum normalized 
   //              correlation between the spectrum of the
   //              envelope of the k-th spectral band of 
   //              x{n} and the spectrum of the envelope 
   //              of the k-th spectral band of x{m}
   //       fo:    Vector containing the center frequencies
   //              of the band filters    
   //       X:     1-D cell array of band-filtered signals
   //              such that X{n}(i,:) is the filtered
   //              version of x{n} through the i-th band
   //              filter. Each signal keeps the original
   //              sample rate
   //       Xenv:  1-D cell array of downsampled versions 
   //              of the envelopes of the band-filtered
   //              signals such that Xenv{n}(i,:) is the
   //              envelope of X{n}(i.:). All envelopes 
   //              have the same sample rate of 200 Hz
   //       tenv:  Envelope time vector in s
   //
   // Signals x{n} are filtered with a filter bank, then 
   // the energy envelopes are computed by squaring the 
   // output of each filter, filtering by a low-pass filter
   // and taking the square root. After downsampling to a
   // sample rate Fsenv = 200 Hz, the spectra of all the 
   // envelopes are computed and the correlations between 
   // the spectra of corresponding envelopes of each pair of
   // signals are computed. The maxima are the components 
   // of Xrho.  
   // 
   // The function is useful to compare signals from the
   // potential semantic content point of view. It turns 
   // out that the correlation between spectra of envelopes
   // is better suited than the correlation between the
   // envelopes themselves.
   // 
   // ----------------------------------
   // Federico Miyara
   // 2006-08-15
   // 2025-09-30
   // 2025-10-06

   // -----------------------------------------------------------------
   // For testing purposes --DELETE--
   // -----------------------------------------------------------------
   if 1==2
      // Duration
      T = 10;
      // Generate white noise normalized to a
      // peak factor of 5 
      Fs(1) = 44100;
      x{1} = rand(1, T*Fs(1), "normal");
      x{1} = x{1} / stdev(x{1}) / 5;
      t{1} = (0:T*Fs(1)-1)/Fs(1);

      // Load speech signal files
      // Sample rates might differ as far as they
      // ensure equal intelligibility conditions 
      path = "s:\balanceadas\"
      file1 = "semántico.wav"
      [x{2}, Fs(2)] = wavread(path + file1);
      x{2} = x{2}(1, 1:T*Fs(2));
      t{2} = (0:T*Fs(2)-1)/Fs(2);

      path = "s:\balanceadas\"
      file2 = 'Memoria.wav' 
      [x{3}, Fs(3)] = wavread(path + file2);
      x{3} = x{3}(1, 1:T*Fs(3));
      t{3} = (0:T*Fs(3)-1)/Fs(3);

      path = "R:\Mis documentos\Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\"
      file3 = "Enigma (FSM).wav"
      [x{4}, Fs(4)] = wavread(path + file3);
      x{4} = x{4}(1, 1:T*Fs(4));
      t{4} = (0:T*Fs(4)-1)/Fs(4);
      // Tipo de bandas: críticas
      typ = 0;    
   end
   // -----------------------------------------------------------------
   
   rhs = argn(2);
   
   if rhs<5
      // Default envelope sample rate
      Fsenv = 200;
   end
   if rhs<4
      // Default duration 
      T = 10;
   end
   if rhs<3
      // Default type of band: critical
      typ = 0;
   end
   
   // Number of signals
   N = length(x);
   
   for n=1:N
      if type(x{n})==10
         // If the cell content is a string then 
         // consider it a file path of a wav file
         // and load it
         [x{n}, Fs(n)] = wavread(x{n});
      end         
      // Reduce to mono
      if size(x{n},1)>1
         x{n} = x{n}(1, 1:T*Fs(n));
      end
   end
   
   // Load function bands()
   exec bands.sci;
   
   // Time constant for envelopes
   tau = 0.01;
   // Do not downsample in the process of band
   // filtering 
   nd = 1;  
   
   // Apply filter bank and get envelopes
   tic
   for n = 1:N
         [X{n}, Xenv{n}, fo] = bands(x{n}, typ, Fs(n), tau, nd);
   end
   toc
   
   // The envelopes will be downsampled to Fsenv to
   // reduce computational load

   // Downsample envelope using intdec()
   tic
   for n=1:N
      for k=1:length(fo)
         XXenv{n}(k,:) = intdec(Xenv{n}(k,:), Fsenv/Fs(n));
      end
   end
   toc
   Xenv = XXenv;
   clear XXenv
   
   
   // We shall correlate the spectra of the envelopes
   // of homologous bands of two signals. These 
   // envelopes have been downsampled at a sample
   // rate of 200 Hz, allowing to include frequencies
   // up to 100 Hz. However, this is still too high 
   // to represent perceptible envelope features so we
   // will set to 25 Hz the maximum relevant frequency
   
   // Maximum relevant frequency
   Fenvmax = 25;
   // FFT frequency bin width
   deltaf = 1/T
   // Maximum relevant spectral line
   hmax = Fenvmax/deltaf;
   // Vector of indices for the discrete frequencies 
   // to keep
   h = 1:hmax;

   // Compute the FFT spectra of the envelopes of
   // all signals, remove their mean values and 
   // correlate pairwise the spectra truncated to 
   // Fenvmax
   tic
   for i = 1:length(fo)
      for n = 1:N
         // Compute the FFT of the i-th band of Xenv{n}
         // with its mean removed
         // NOTE: The mean of the envelope is much larger
         //       than the other spectral components of
         //       the envelope. It is not relevant but 
         //       could affect the correlation 
         A{n,i} = abs(fft(Xenv{n}(i,:) - mean(Xenv{n}(i,:))));
         
         // Remove the mean value of each spectrum to 
         // increase sensitivity of the correlation test
         B{n,i} = A{n,i} - mean(A{n,i});
      end
      
      // Correlation between pairs of spectra of 
      // spectral envelopes of signals
      for n=1:N
         for m = n:N
            C = xcorr(B{n,i}(h), B{m,i}(h));
            C = C / sqrt(sum(B{n,i}(h).^2) * sum(B{m,i}(h).^2));
            Xrho(n,m,i) = max(C);
            Xrho(m,n,i) = Xrho(n,m,i);
         end 
      end
      
      // Downsampled envelopes time vector
      tenv = (0:size(Xenv{1},2)-1)/Fsenv;

   end
   toc

   // disp([round(fo(1:20)) ..
   //       rho(1,2,:)(:) rho(1,3,:)(:) rho(1,4,:)(:) ..
   //       rho(2,3,:)(:) rho(2,4,:)(:) rho(3,4,:)(:)])

endfunction



