function y = notched_noise(f, BW, T, Fs, typ, N)
   // Notched noise generator
   // 
   // Generates a white noise from which a band of 
   // given bandwidth around a given frequecy has 
   // been filtered out with a band-reject filter
   // with steep edges
   // 
   // Usage:
   //       y = notched_noise(f, BW, T, Fs, typ)
   // where
   //       f:   Centre frequency in Hz
   //       BW:  Bandwidth in Hz
   //       T:   Duration in s
   //       Fs:  Sample rate in Hz
   //       typ: Type of noise; if "p", uses perfect
   //            white noise, otherwise, uses random
   //            noise
   //       N:   Window size       
   //
   // NOTE: The perfect noise version generates a very
   //       deep spectral well, at the expense of being
   //       noticeably periodic. Periodicity is hardly
   //       perceived for N >= 65536. The random noise
   //       version is almost truly random, but the 
   //       spectral well is not so deep due to spectral
   //       leakage of the overlap-add technique in the 
   //       FFT filter algorithm   
   //
   // Needs critband(), whitenoise_perfect(), whitenoise()
   //       and fftfilter()
   //
   // ----------------------------------
   // Author: Federico Miyara
   // Date:   2024-09-29 

   rhs = argn(2);
   
   // Missing arguments
   if rhs<6
      // Filtering window size
      N = 32768;   
   end
   if rhs<5
      typ = "random";
   end   
   if rhs<4 
      Fs = 44100;
   end
   if rhs<3
      T = 5
   end
   if rhs<2
      // If no bandwidth is provided, the tritical
      // bandwidth is used 
      exec critband.sci;
      BW = critband(f, 0)
   end
   

 
   // Filtering window
   w = ones(1,N/2);
   fi = (sqrt(BW^2 + 4*f^2) - BW)/2
   fs = (sqrt(BW^2 + 4*f^2) + BW)/2
   ki = ceil(fi*N/Fs) + 1
   ks = floor(fs*N/Fs) + 1
   w(ki:ks) = 0;
   W = [w, 0, conj(w($:-1:2))];
   
   if typ=="p" | typ=="perfect"
      // Generate one window of perfect white noise
      exec whitenoise_perfect.sci;
      xo = whitenoise_perfect(N/Fs, Fs, N);
      // Apply filter in the frequency domain to
      // a single window and replicate
      yo = ifft(fft(xo).*W);
      M = T*Fs;
      m = floor(M/N);
      y = repmat(yo, 1, m);
      y = [y, yo(1:M-m*N)];
   else
      // Generate white noise
      exec whitenoise.sci;
      x = whitenoise(T, Fs);
      // Apply FFT filter
      exec fftfilter.sci;
      y = fftfilter(x, w, N);
   end
   
endfunction
