function y = psicoenhancer(x, f, L, e, Fs)
   // Sound signal psichoacoustic enhancer
   // This effect takes a sound signal, applies a nonlinear
   // distortion to the spectral components above frequency
   // f and then adds this signal to the dry, undistorted
   // full-band signal.
   //
   // Usage:
   //        y = psicoenhancer(x, f, L, e, Fs)
   //
   // where
   //        x:  original signal
   //        f:  cutoff frequency of highpass filter in Hz
   //        L:  mix level of wet signal (dry = 0 dB)
   //        e:  odd (1) or even (2) harmonics
   //        Fs: sampling rate in Hz
   //
   // f, L, e and Fs are optional, being their default values 
   // 1000, -20, 2 and 44100.
   //
   // The function distorts the input signal so that high
   // frequency harmonics are added. In contrast to a plain
   // distortion, which would add also noticeable and 
   // disagreable low frequency distortion products, the
   // enhancer adds spectral components contributing to 
   // brightness.
   //  
   // -------------------------------
   // Author: Federico Miyara
   // Date:   2010-06-20
   //         2012-08-23
   //         2015-05-31
   //         2022-08-20

   // Default values
   rhs = argn(2)
   if rhs < 5
      Fs = 44100;
   end
   if rhs < 4
      e = 2;
   end
   if rhs < 3
      L = -20;
   end
   if rhs < 2
      f = 1000;
   end

   // High-pass filtering window
   N = 4096;
   wh = zeros(1, N/2);
   wh(round(f/Fs*N):$) = 1; 

   // High-pass filtering
   exec fftfilter.sci;
   xh = fftfilter(x, wh, N);

   // Normalization of high-pass signal
   xh = xh/(max(abs(xh(100:$))));

   // Distortion of high-pass signal (wet signal)
   if e == 1
      // Odd-order harmonics
      xhd = sign(xh).*(1 - exp(-10*abs(xh)));
   else
      // Even-order harmonics
      xhd = (1 - exp(-5*abs(xh))).^2;
   end

   // Remove any low frequency distortion products
   xhd = fftfilter(xhd, wh, N);

   // Addition of the dry (original) and wet signals
   y = x + 10^(L/20) * xhd;

   //figure(1)
   //plot(t, x, t, xh, t, xhd, t, y)

   //f = [0:N-1]*Fs/N;
   //figure(2)
   //plot(f, abs(fft(x,N)), f, abs(fft(y,N)))

endfunction
