// ====================================================================================================== //
// Author      :: Sumit Adhikari
// Contact     :: adhikari@ict.tuwien.ac.at,adhikari@ieee.org
// Affiliation :: ICT, Vienna University of Technology,
//             :: Gusshausstrasse 27-29/E384, 1040 Vienna
// ====================================================================================================== //

// ------------------------------------ Initial Setup --------------------------------------------------- //
   clear all          ;
   funcprot(0)        ;
   stacksize('max')   ;
   show_issue     = 0 ;

   function [ACC] = power(B)
      ACC = 0          ;
      for i=1:size(B,"*")
         if sign(B(i)) == -1 then bitsign = 1;
         else                     bitsign = 0;
         end
         ACC = ACC + ceil(log2(B(i))) + bitsign ;
      end
   endfunction

   function [] = getArray(B)
      for i=1:size(B,"*")
         printf(" B(%02d) = %0d\n",i,B(i));
      end
   endfunction

   function [] = getArraySorted(B)
      B = mtlb_sort(B);
      for i=1:size(B,"*")
         printf(" B(%02d) = %0d\n",i,B(i));
      end
   endfunction

   
   function [A] = invert(B)
      N = size(B,"*");
      for i=1:N
         A(i) = B(N-i+1);
       end
   endfunction 


   BW  = 1024      ;
   Fc  = 120.0     ;
   
// --------------------------------- Filter --------------------------------- //
   M2        = 11 ;
   Hiir      = iir(M2,'lp','ellip',[Fc/BW 150/BW],[.01 .0001]);
   [hiir,fr] = frmag(Hiir,BW);
   clf(2);xset('window',2);plot(BW*fr',20*log10(hiir),"red");

// ------------------------------------ Filter Quantization ----------------------------------------- //
   Nbits  = 28 ;
   W      = (2^Nbits - 1) ;

   b      = fix(round(W*coeff(numer(Hiir))));
   a      = fix(round(W*coeff(denom(Hiir))));
   
   Hn     = fix(round(W*numer(Hiir)));
   Hd     = fix(round(W*denom(Hiir)));


   Hiir_q = Hn/Hd ;


   b      = coeff(numer(Hiir_q));
   a      = coeff(denom(Hiir_q));


// ============ Show-Up Issue ===================== //
if show_issue == 1 then

   z = poly(0,'z');
   Hnumer = syslin('c',(157237*z^11 -573356*z^10 +1283737*z^9 -1748438*z^8 +1626847*z^7 -622496*z^6 -622496*z^5 +1626847*z^4 -1748438*z^3 +1283737*z^2 -573356*z^1 + 157237 ),(z^16)); 
   Hdenom = syslin('c',(268435455*z^11 +2118480371*z^10 -173903078*z^9 +1115050201*z^8 -248360047*z^7 +1491494717*z^6 -1417205199*z^5 +1076506307*z^4 -1042913473*z^3 +383618516*z^2 +801201069*z^1 -77190481 ),(z^16));

   Hiir_q2 = Hnumer/Hdenom ;
   b      = coeff(numer(Hiir_q2));
   a      = coeff(denom(Hiir_q2));

   Hiir_q =  Hiir_q2

end


   [hiir_q,fr] = frmag(Hiir_q,BW);
   plot(BW*fr',20*log10(hiir_q),"green");
   xtitle('Response - Original And Quantized','Frequency (Hz)','Gain (dB)','Z axis') ;mtlb_grid("on");

// --------------------------------- Pole-Zero Plot for Hiir ------------------------------------------ //
   clf(4);xset('window',4);plzr(Hiir)  ; xtitle('Pole-Zero Plot for Hiir')  ; mtlb_grid("on");
   clf(5);xset('window',5);plzr(Hiir_q); xtitle('Pole-Zero Plot for Quantized Hiir'); mtlb_grid("on");