% sequential direct form II IIR implementation test
% 2 delays
% 6 multiplications (worst case)
clear all
g = 0.0421;
b = [1 -0.58570773849879676 1].*g;
a = [1 -1.6888219088693734 0.75566935461656326];

input = zeros(1000,1);
input(1) = 1;
output = zeros(1000,1);
f = linspace(-1,1,length(input)); % normalized frequency

int = 0;
intd = 0;
intdd = 0;

for n=1:length(input)
   intdd = intd;
   intd = int;
   int = a(1)*input(n) - a(2)*intd - a(3)*intdd;
   output(n) = b(1)*int + b(2)*intd + b(3)*intdd;
end

subplot(3,1,1)
stem(input)
hold on
grid on
stem(output,'r')
xlim([0 30])
title('Impulse Response')
xlabel('sample no.')
ylabel('amplitude')

subplot(3,1,2)
plot(f, 20*fftshift(log10(abs(fft(input)))))
hold on
grid on
plot(f, 20*fftshift(log10(abs(fft(output)))),'r')
title('Frequency Plot');
xlabel('normalized frequency')
ylabel('magnitude in dB')

subplot(3,1,3)
plot(f, fftshift(angle(fft(input))))
hold on
grid on
plot(f, fftshift(angle(fft(output))),'r')
title('Phase Plot');
xlabel('normalized frequency');
ylabel('phase in rad')

