Three times is a charm. Maybe without the raw tags this time:
clc
clear all
pkg load signal
pkg load control
graphics_toolkit fltk
a0 = 1.096778649043814
a1 = -0.2986467797702343
a2 = 0.43403424153896397
b1 = -0.2986467797702343
b2 = 0.5308128905827781
zeros3=[a0 a1 a2]
poles3=[1.0 b1 b2]
a0 = 0.24807653063671162
a1 = 0.21331109909099844
a2 = 0.19896775914001172
b1 = 0.21331109909099844
b2 = -0.5529557102232766
zeros30=[a0 a1 a2]
poles30=[1.0 b1 b2]
figure
freqz(conv(zeros3,zeros30),conv(poles3,poles30),1024,48000.0)
sysd=tf(zeros3,poles3,1/48000.0)
combi_coeff_zeros=conv(zeros3,zeros30);
combi_coeff_poles=conv(poles3,poles30);
function height=calch(radius,combi_coeff_zeros, combi_coeff_poles)
angl=[0:pi/64:pi];
height=repelem(0,length(angl))
% 64 steps n 1..64 to go round the circle
for n=1:(length(angl))
re=cos(angl(n));
im=sin(angl(n));
z=abs(polyval(combi_coeff_zeros,radius*complex(re,im)));
p=abs(polyval(combi_coeff_poles,radius*complex(re,im)));
height(n) = 20*log10(z/p);
end
end
steps=100;
stepshalf=idivide(steps,2);
%height=zeros(steps,steps);
height=zeros(steps,stepshalf);
x=1;
rangeunit=([0.5-stepshalf:stepshalf-0.5]/stepshalf);
rangeunithalf=([0.0:stepshalf-0.5]/stepshalf);
for re=rangeunit
y=1;
for im=rangeunithalf
z=abs(polyval(combi_coeff_zeros,complex(re,im)));
p=abs(polyval(combi_coeff_poles,complex(re,im)));
height(x,y) = 20*log10(z/p);
y++;
end
x++;
end
figure
mesh(rangeunithalf,rangeunit,height);
hold on
%pole-zero plot
ze=roots(combi_coeff_zeros);
po=roots(combi_coeff_poles);
vert=[-40:5:40];
angl=[0:pi/128:pi];
for v=vert
plot3(abs(imag(ze)),real(ze),repelem(v,length(ze)),'ok')
plot3(abs(imag(po)),real(po),repelem(v,length(po)),'xr')
%draw circle:
%plot3(sin(angl),cos(angl),repelem(v,length(angl)),'r')
end
x=[pi:-pi/32:-pi]/pi;
an=[pi/2:-pi/64:-pi/2];
h=repelem(0,length(x));
h=calch(1.0,combi_coeff_zeros,combi_coeff_poles);
% the frequency response of the filter in blue, projected outside the circle
plot3(repelem(1,length(x)),x,h);
% the intersection with the circle, in black, it is the frequency response
wrapped around the tube!
plot3(cos(an),sin(an),h,'k');
--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Thanks