function y = pinknoise(T, Fs, N, mtd)
   // Pink noise generator
   //
   // Usage:
   //       y = pinknoise(T, Fs, N, mtd)
   //
   // where T:   duration in s
   //       Fs:  sampling rate in Hz
   //       N:   FFT filtering window size in samples
   //       mtd: method
   //            1: the FFT spectrum is 1/sqrt(n)
   //            2: the spectral density is 1/f
   //       y:   pink noise output row vector 
   //
   // mtd, N and Fs are optional. The default values are,
   // respectively, 1, 4096 and 44100. N must be an even 
   // number but it is strongly recommended that it be a
   // power of 2. For N = 16384 the spectral slope has an
   // error smaller than 1 dB in the whole audio band,
   // 
   // The function works by FFT-filtering a white noise
   // with a filtering window with a 3 dB/oct decay. 
   // The root mean square value is normalised to 1. The
   // -3 dB/oct is kept steadily up to the Nyquist limit.
   // Needs fftfilter.sci in the same folder.
   //
   // Se also: whitenoise, browniannoise, bandnoise
   //
   // ------------------------------------
   // Author: Federico Miyara
   // Date:   2010-06-07
   //         2019-05-24
   //         2020-08-29
   //         2023-03-07

   // Provide default values for missing arguments
   rhs = argn(2);
   if rhs<4
      mtd = 1;
   end
   if rhs<3
      N = 4096;
   end
   if rhs < 2
      Fs = 44100;
   end

   // Filtering window design
   // In order to generate pink noise from white noise
   // the filtering window must have a slope of -3 dB/oct.
   // This is accomplished by letting W be inversely 
   // proportional to the square root of frequency
   n = 1:N/2;
   mtd = 1
   select mtd
   case 1
      // Approximate filtering window
      // Yields a 1/sqrt(n) FFT slope, which is 
      /// not exactly equivalent to a 1/f 
      // spectral density, especially at the
      // low frequency end
      w = 1 ./ sqrt(n);
   case 2
      // Slightly more exact filtering window
      // Yields 1/f spectral density
      w = sqrt(2 * log((n+1)./n) / log(N/2));
   end

   // White noise gaussian noise as column vector 
   y = grand(1, 1 + floor(Fs*T), "nor", 0, 1);


   // Load function
   flag = exists("fftfilter")
   if ~flag
      exec fftfilter.sci;
   end

   // Signal filtering
   y = fftfilter(y, w, N);

   // Unload function if previously not loaded
   if ~flag
      clear fftfilter;
   end

   // Remove DC component
   y = y - mean(y);

   // Normalization to get RMS value equal to 1.
   y = y / sqrt(mean(y.^2));

endfunction
