1 % by frisch 
 2 
 3 clear;
 4 n=1024; 
 5 f=linspace(0,1,n);
 6 w=pi*f; 
 7  
 8 % B = [1.0770389918072; -1.03115724017531; 1.261079950480885];
 9 % A = [1; -1.98606828810359; 0.992414110487639];
10 B = [0.03; 0.6605; 0.121; 0.0605; 0.03]; 
11 A = [1; -1.194; 0.436]; 
12 % compute frequency response of Hd(w): 
13 z=exp(i*w);
14 z_numerator =[ones(1,n); z.^(-1); z.^(-2);z.^(-3);z.^(-4);]; 
15 z_denominator =[ones(1,n); z.^(-1); z.^(-2)]; 
16 Hdw=(B'*z_numerator)./(A'*z_denominator); 
17 % plot Group Delay:
18 dpHdw = diff(angle(Hdw))./diff(w); 
19 plot(f(2:end),-dpHdw,'r') 
20 title('Group Delay','FontName','Times','FontSize',12,'FontWeight','bold'); 
21 xlabel('Frequency [Hz]','FontName','Times') 
22 ylabel('Group Delay [s]','FontName','Times') 
23 %axis([0 1 1 4])  
24 grid on