Re: 3d pole-zero plot in z-domain

[address@hidden]

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