function [y, t] = imp_resp(H, Fs)
   // Compute impulse response from frequency response
   // 
   // Usage:
   //       [y, t] = imp_resp(H, Fs)
   // where
   //       H:  complex frequency response defined at
   //           a number of equally spaced frequencies,
   //           the first one is 0 (DC) and the last
   //           one is Fs/2*(length(h) - 1)/length(h) 
   //       Fs: sample rate in Hz
   //       y:  impulse response vector
   //       t:  time vector in s
   // 
   // NOTES: 1) Usually H is computed from a transfer
   //           function, but could also be the result
   //           of measurement of a real system
   //        2) Fs is optional, its default value is
   //           44100
   //        3) The length of the impulse response is
   //           limited to the spectral data available,
   //           i.e., 2*length(H). To get a longer
   //           impulse response either provide more
   //           spectral data or else interpolate the
   //           available data
   //     
   // ------------------------
   // Author: Federico Miyara 
   // Date:   2023-04-13
   //         2024-07-23
   
   // test data --delete--
   if 1==2
      N = 65536;
      Fs = 44100;
      f = (0:N/2-1)*Fs/N;
      fo = 1000;
      H = 1./(1 + %i*f/fo);
   end
   
   // Set default value for missing argument
   if argn(2)<2
      Fs = 44100;
   end
   
   // The transfer function is equal to the spectrum 
   // of the impulse response. Extend the frequency
   // response to frequencies above Nyquist taking
   // into account that for real signals there is a
   // conjugate symmetry condition for the FFT.  
   H = [H, 0, conj(H($:-1:2))];
   
   // Get the inverse fast Fourier Transform to
   // retrieve the unit pulse response  
   y = ifft(H);
   
   // Get the impulse response denormalizing by Fs
   // since the unit pulse is the discrete equivalent
   // to the impulse response divided by Fs
   y = Fs*y;
   
   // Time vector
   t = (0:length(H)-1)/Fs;     
        
endfunction
