function x = bandnoise(f1, f2, T, Fs, Xef, N, method)
   // This function generates band limited random noise
   //
   // Usage:
   //       x = bandnoise(f1, f2, T, Fs, Xef, N, method)
   // 
   // where
   //       f1:     low frequency cutoff in Hz
   //       f2:     high frequency cutoff in Hz
   //       T:      duration in s
   //       Fs:     sampling rate in Hz
   //       Xef:    effective value of noise
   //       N:      filtering window size (power of 2)
   //       x:      row vector containing noise samples
   //       method: "rand": derived from random white noise
   //               "perf": spectrally-perfect band noise    
   //
   // T, Fs, Xef and N are optional, being their default 
   // values 1, 44100, 1 and 4096 respectively.
   //
   // While N may be any even number, the algorithm is faster 
   // when it is a power of 2. If f1 = f2, a single spectral
   // line is kept. The same is true if f2 - f1 << Fs/N; 
   //
   // The function generates a white noise and then filters it 
   // with a filtering window that is 1 at the in-band
   // frequencies and 0 elsewhere. Filtering is performed with 
   // a fast FFT-based overlap-add algorithm.
   //
   // NOTE: The method can be "rand" (default), i.e. the band
   //       noise is obtained filtering a long-term average 
   //       white noise, or "perf", filtering a spectrally
   //       perfect white noise. This means that when analyzed
   //       with an FFT of size N it is perfectly flat in 
   //       magnitude but random in phase. If choosing "perf"
   //       N should be 65536 or greater to avoid an audible
   //       repeated pattern         
   //
   // ----------------------------
   // Author: Federico Miyara
   // Date:   2010-06-03
   //         2020-01-12
   //         2020-08-29
   //         2023-03-07
   //         2024-11-03

   // Number of arguments
   rhs = argn(2)

   if rhs<7
      method = "rand";
   end
   if rhs<6
      N = 4096;
   end
   if rhs<5
      Xef = 1;
   end
   if rhs<4
      Fs = 44100;
   end
   if rhs<3
      T = 1;
   end

   // White noise generation
   if method=="perf"
      exec whitenoise_perfect.sci;
      x = whitenoise_perfect(T, Fs, N);
   else
      x = grand(1, floor(T*Fs) + 1, "nor", 0, 1);
   end

   // Design of filtering window
   // Window initialisation
   W = zeros(1, N/2);
   // Index of low frequency cutoff
   k1 = max(1, round(f1/Fs*N));
   // Index of high frequency cutoff
   k2 = round(f2/Fs*N);
   // Turn on in-band window values
   W(k1:k2) = ones(1, k2-k1+1);

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

   // Filter to get band noise
   x = fftfilter(x, W, N);

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

   // Remove DCcomponent
   x = x - mean(x);

   // Normalisation to comply with prescribed Xef
   x = x / sqrt(mean(x.^2)) * Xef;

endfunction
