function [z, t, f] = spectrogram(x, N, M, Fs, w, grph, cmap, DR)
   // Spectrogram of a signal
   //
   // Usage
   //       [z, t, f] = spectrogram(x, N, M, Fs, w, grph, cmap, DR)
   // where
   //       x:    signal vector
   //       N:    FFT window length. Default: 1024
   //       M:    offset between FFT windows in samples
   //             Default: 512 
   //       Fs:   sample rate in Hz. Default: 44100
   //       w:    windowing function; accepted values are
   //             "boxcar", "hann", "hamming", "blackman", 
   //             "bharris", "flattop" and a custom vector
   //             containing N components of the window.
   //             Default: "hann"
   //       grph: graphic output; accepted values are
   //             "lin", "lincmap", "log", "logcmap",
   //             "newlin", "newlincmap", "newlog",
   //             "newlogcmap", "none"; lin/log refers 
   //             to linear or logarithmic distribution
   //             of magnitude along the color map, cmap
   //             to the inclusion or not of a color map, 
   //             and new to the creation of a new 
   //             window. The default is "newlog".
   //       cmap: a string indicating any of the accepted
   //             color maps (such as "jet") or any matrix 
   //             with 3 columns containing the RGB
   //             description of the colors. Default: "jet"
   //             (see Predefined Colormaps)
   //       DR:   dynamic range in dB, max 0 dB min -DR dB. 
   //             Default: 100 dB
   //       z:    matrix with N/2 rows and
   //             floor((length(x) - N)/M + 1) columns, 
   //             where each column contains the absolute
   //             value of the spectrum beginning at each 
   //             multiple of M
   //       t:    vector containing center times in s of
   //             each FFT window
   //       f:    vector containing the values of the
   //             FFT frequency bins in Hz
   //
   // The numerical data provided as output arguments can 
   // be used to get information about spectral features 
   // or to obtain other graphic presentations, such as 
   // time slices (selecting specific times) or 3D, surface
   // or waterfall representations.
   // z is scaled so that for a sine wave whose frequency
   // is coincident with the k-th spectral bin z(k) is
   // equal to its amplitude.  
   // NOTE: The function is originally intended to get
   //       spectrograms of sound signals whose maximum 
   //       amplitude is 1. This condition can be relaxed
   //       if it can be guaranteed that the amplitudes
   //       of all the spectral components is smaller than 
   //       or equal to 1. Otherwise, the values exceeding
   //       1 will be represented with the same colour as
   //       1. This concern applies only to the graphic 
   //       output.  
   //
   // -----------------------------------------------
   // Author: Federico Miyara
   // Date:   2021-07-10
   //         2023-10-29
   //         2024-05-28
   //         2025-01-19
   //         2025-02-01
   //         2025-03-13
   //         2025-03-23

   // Number of colors
   C = 64;

   // Argument handling
   rhs = argn(2);
   lhs = argn(1);

   // NOTE: The alias _ of gettext has been used for
   //       error messages  
   if lhs>3
      // More than 3 output parameters
      msg = _("%s: Wrong number of output arguments: %d to %d expected.\n")
      error(msprintf(msg, "spectrogram", 0,3))
   end
   if rhs==0 | rhs>8
      msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
      error(msprintf(msg, "spectrogram", 1,8))
   end

   if rhs<8
      // Unspecified dynamic range
      DR = 100;
   else
      // Check dynamic range
      if type(DR)<>1
         msg = _("%s: Argument #%d: Decimal number expected.\n")
         error(msprintf(msg, "spectrogram", 8))
      elseif DR<=0
         msg = _("%s: Argument #%d: Must be > %d.\n")
         error(msprintf(msg, "spectrogram", 8, 0))
      end
   end

   if rhs<7
      // Unspecified colormap
      cmap = jet(C);
   else
      // Check colormap
      errmsg = "Unsupported color map. Using jet(" + string(C) + ")";
      if type(cmap)==10
         // cmap is a string with predefined colormap.
         // Accepted color maps are
         cmaps = ["autumn" "bone" "cool" "copper" "gray" "hot" ..
         "hsv" "jet" "ocean" "parula" "pink" "rainbow" ..
         "spring" "summer" "white" "winter" "blues" ..
         "BrBG" "BuGn" "BuPu" "coolwarm" "GnBu" ..
         "greens" "greys" "oranges" "OrRd" "PiYG" ..
         "PRGn" "PuBu" "PuBuGn" "PuOr" "PuRd" ..
         "purples" "RdPu" "RdYlBu" "RdYlGn" "reds" ..
         "spectral" "RdBu" "RdGy" "YlGn" "YlGnBu" ..
         "YlOrBr" "YlOrRd"];
         if or(cmap==cmaps)
            // Get cmap from accepted color maps
            execstr("cmap = " + cmap + "(C)");
         else
            // cmap is an unknown string
            warning(errmsg);
            cmap = jet(C);
         end
      elseif type(cmap)==1
         // cmap is an array of reals
         if isreal(cmap) 
            if size(cmap)(2)==3 & max(cmap)<=1 & min(cmap)>=0
               // cmap has already the final form. Do nothing
            else
               warning(errmsg);
               cmap = jet(C);
            end
         else
            // cmap is complex
            warning(errmsg);
            cmap = jet(C);
         end
      else
         warning(errmsg);
         cmap = jet(C);   
      end
   end

   if rhs<6
      // Unspecified graphic output type
      grph = "log";
   else
      // Check graphic output type
      if type(grph)==10
         // grph is a string
         grphs = ["lin" "lincmap" "log" "logcmap" "newlin" ..
         "newlincmap" "newlog" "newlogcmap" "none"];
         if ~or(grph==grphs)
            // grph is an unknown string
            msg = _("%s: Argument #%d: Unsupported graphic mode.\n")
            error(msprintf(msg, "spectrogram", 6));
         end
      else
         msg = _("%s: Argument #%d: Unsupported graphic mode.\n")
         error(msprintf(msg, "spectrogram", 6));
      end
   end

   if rhs<5
      // Unspecified window type
      w = "hann";
   else
      // Check window type
      if type(w)==10
         winds = ["boxcar" "hann" "hamming" "blackman" ...
         "bharris" "flattop"];
         if ~or(w==winds)
            msg = _("%s: Argument #%d: Unsupported window.\n")
            error(msprintf(msg, "spectrogram", 5));
         end
      elseif type(w)==1
         if and(size(w)==[1, N]) | and(size(w)==[N, 1])
            // w is the right size, do nothing 
         else
            // w has incompatible dimensions
            msg = _("%s: Argument #%d: Incorrect window size.\n")
            error(msprintf(msg, "spectrogram", 5));
         end
      else
         msg = _("%s: Argument #%d: Unsupported window.\n")
         error(msprintf(msg, "spectrogram", 5));
      end

   end

   if rhs<4
      // Unspecified sample rate
      Fs = 44100;
   else
      // Check sample rate
      if type(Fs)<>1
         msg = _("%s: Argument #%d: Decimal number expected.\n")
         error(msprintf(msg, "spectrogram", 4))
      end
   end

   if rhs<3
      // Unspecified offset between windows
      M = 512;
   else
      // Check offset between windows
      if type(M)<>1
         msg = _("%s: Argument #%d: Decimal number expected.\n")
         error(msprintf(msg, "spectrogram", 3))
      end      
   end

   if rhs<2
      // Unspecified window size
      N = 1024;
   else
      // Check window size
      if type(N)<>1
         msg = _("%s: Argument #%d: Decimal number expected.\n")
         error(msprintf(msg, "spectrogram", 2))
      end      
   end

   // Check input signal
   if type(x)<>1
      msg = _("%s: Argument #%d: Decimal number expected.\n")
      error(msprintf(msg, "spectrogram", 1))
   end



   // Data for testing purposes --DELETE--
   if 1==2
      N = 1024
      M = 256
      Fs = 44100
      w = "hann"
      grph = "logcmap"
      cmap = "jet"
      tt = 0:1/Fs:2;
      fo = 8000
      fm = 0.5
      df = 4000
      ff = fo + df*sin(2*%pi*fm*tt);
      phi = 2*%pi*cumsum(f)/Fs;
      caso = 2
      select caso
      case 1
         // Two harmonics with aliasing
         x = 0.4*sin(phi) + 0.3*sin(2*phi);
      case 2
         // Two harmonics with aliasing
         // and distortion
         x = sin(phi) + sin(2*phi);  
         x = min(x, ones(x));          
      end
      wavwrite(x,"d:\hello\fmod.wav");
   end



   // Ensure the signal is a column vector
   x = x(:);

   // Original signal length
   L = length(x);
   // Number of FFT windows
   q = floor((L+N)/M);
   // Signal length after zero-padding
   LL = (q+1)*M + N;
   // Add zero-padding at the beginning and end
   xx = [zeros(N,1); x; zeros(LL-L-N, 1)]; 
   //xx = [zeros(N/2,1); x; zeros(N/2 - (L - floor(2*L/N)*N/2), 1)]; 

   // Number of FFT windows 
   Q = floor((length(xx) - N)/M);
   //Q = floor((length(x) - N)/M + 1);

   // Initialize z
   z = zeros(N/2, Q);

   // Generate a matrix whose columns are
   // segments with length N of the signal
   // starting at sample indices multiple 
   // of M
   y = zeros(N, Q);
   for k=1:Q
      y(:,k) = xx((k-1)*M+1:(k-1)*M+N);
   end

   // Compute window function
   if type(w)==1
      w = w(:);
   else
      select w
      case "boxcar"
         w = ones(1,N)';
      case "hann"
         w = 0.5 - 0.5*cos(2*%pi*[0:N-1]/N)';
      case "hamming"
         w = 0.54 - 0.46*cos(2*%pi*[0:N-1]/N)';
      case "blackman"
         w = 0.42 - 0.5*cos(2*%pi*[0:N-1]/N)' + ...
         0.08*cos(4*%pi*[0:N-1]/N)'; 
      case "bharris"
         w = 0.35875 - 0.48829*cos(2*%pi*[0:N-1]/N)' + ... 
         0.14128*cos(4*%pi*[0:N-1]/N)' - ...
         0.01168*cos(6*%pi*[0:N-1]/N)';
      case "flattop"
         w = 1 - 1.90796*cos(2*%pi*[0:N-1]/N)' + ...
         1.07349*cos(4*%pi*[0:N-1]/N)' - ...
         0.18199*cos(6*%pi*[0:N-1]/N)';
      end
   end
   // Replicate window
   ww = repmat(w, 1, Q);


   // FFT of the windowed signal along columns
   // and conversion to unilateral spectrum
   z = abs(fft(y.*ww, -1, 1))(1:N/2,:)*2/N;
   z(1,:) = z(1,:)/2;

   // Window attenuation compensation
   W = sqrt(mean(w.^2));
   z(2:$,:) = z(2:$,:)/W;
   z(1,:) = z(1,:)/mean(w);

   // Limit dynamic range
   // NOTE: This is necessary to avoid, in log
   //       scale, problems due to zero values
   z = max(z, max(z)*10^(-DR/20));

   // Starting times of the FFT windows 
   t = ((0:Q-1)*M - N/2)/Fs;

   // Frequencies of the FFT bins
   f = (0:N/2-1)*Fs/N;

   // Remove spectra outside signal duration
   T = L/Fs;
   I = find((t>=0)&(t<=T));
   t = t(I);
   z = z(:,I);

   // Graphic output
   select grph
   case "lin"
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, z', zminmax=[10^(-DR/20) 1], colout=[-1 -1]);
   case "newlin"
      scf();
      clf();
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, z', zminmax=[10^(-DR/20) 1], colout=[-1 -1]);
   case "lincmap"
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, z', zminmax=[10^(-DR/20) 1], colout=[-1 -1]);
      colorbar(10^(-DR/20), 1)   
   case "newlincmap"
      scf();
      clf();
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, z', zminmax=[10^(-DR/20) 1], colout=[-1 -1]);
      colorbar(10^(-DR/20), 1)   
   case "log"
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
   case "newlog"
      scf()
      clf();
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
   case "logcmap"
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
      colorbar(-DR, 0)    
   case "newlogcmap"
      scf()
      clf();
      gcf().color_map = cmap;
      gcf().background = -2;
      Sgrayplot(t, f, 20*log10(z'), zminmax=[-DR 0], colout=[-1 -1]);
      colorbar(-DR, 0) 
   case "none"  
   end

endfunction

