dfdpzp3.m

很多matlab的源代码

function [num,den] = dfdpzp3(tol)
% DFDPZP3 Interactive digital filter design by pole zero placement.
%
%         [NUM,DEN] = DFDPZP3(TOL) invokes an interactive digital filter
%         design by pole-zero placement with tolerance, TOL, for snapping. 
%         [Default: TOL = 0.02].  Upon completion, NUM and DEN contain the
%         coefficients in descending powers of z.
%
%                   H(z) = NUM(z)/DEN(z)
%
%         Left mouse button places poles, Right mouse button places zeros.
%         <CR> Enter key plots or ends program.
%         Defaults: If no zeros present, NUM = 1. 
%                   If no poles present, DEN = 1.
%
%         Note: THIS ROUTINE IS ONLY FOR MATLAB VERSION 3.5


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if exist('version')==5,disp('This routine is ONLY for MATLAB v3'),return,end
clc,disp('')
disp('Interactive digital filter design by pole zero placement')
disp('Place poles and zeros on the screen using mouse buttons')
disp('L button places poles, R places zeros, Enter key plots or ends prog')
disp('DEFAULTS: If no zeros present, Num=1. If no poles present, Den=1.')
disp(' '),disp('Strike a key to begin'),pause

if nargin==0,tol=0.02;end
j=sqrt(-1);subplot,hold off,polar(2*pi,1),grid,hold on
title('LB:POLES  RB:ZEROS  ENTER KEY:PLOTS')
z=[];p=[];bb=0;while bb ~= 13
[xx,yy,bb]=ginput(1);
if abs(yy)<tol,yy=0;xx=0;end  %%SNAP TO ORIGIN
if abs(yy*sin(xx))<tol,xx=0;end %%SNAP TO REAL VALUE
if xx ~= 0,xx=[xx;-xx];yy=[yy;yy];end

if bb==1,polar(xx,yy,'xc5'),pp=yy.*exp(j*xx);p=[p;pp];end
if bb==2,if abs(1-abs(yy(1)))<tol,yy=0*yy+1;end,polar(xx,yy,'oc12')
zz=yy.*exp(j*xx);z=[z;zz];end
end,

if isempty(z),num=1;else,num=real(poly(z.'));end              %%NEW
if isempty(p),den=1;else,den=real(poly(p.'));end              %%NEW
n=0:100;n0=0*n';n0(1)=1;y=filter(num,den,n0);
w=0:pi/200:pi;zzz=exp(sqrt(-1)*w);h=polyval(num,zzz)./polyval(den,zzz);
hold off,axis('normal');subplot(211)
dtplot(n,y,'.'),grid,title('MENU: LB: DEL  RB: ADD  ENTER KEY: END PROG'),...
xlabel('Impulse response n=0:100'),...
subplot(212),plot(w/2/pi,abs(h)),grid,xlabel('Frequency Response vs F'),...
subplot,hold off
[xxx,yyy,bbb]=ginput(1);

while bbb ~= 13
if bbb==1 | bbb==2,subplot,hold off,polar(2*pi,1),grid,hold on
polar(angle(p),abs(p),'xc5'),polar(angle(z),abs(z),'oc12'),end

if bbb==1,title('DELETE: LB:POLES  RB:ZEROS  ENTER KEY: PLOTS');
bb=1;while bb ==1 | bb==2
[xx,yy,bb]=ginput(1);if bb==1,if ~isempty(p),
%pa=abs(xx-angle(p));%%REPLACED BY BETTER COMPARISON
pa=abs(yy.*exp(j*xx)-p);
k=[];k1=[];i=find(pa==min(pa));
if ~isempty(i),i1=i(1);  %%NEW
if imag(p(i1))~=0,k=find(p==conj(p(i1)));
if ~isempty(k),k1=k(1);end,end,end %NEW
l=[i1 k1];
if (isempty(k1) & length(i)==1)|length(k)==1 %%NNN
polar(angle(p(l)),abs(p(l)),'xi'),
end %NNN
p(l)=[];end,end

if bb==2,if ~isempty(z),
%za=abs(xx-angle(z));%%REPLACED BY BETTER COMPARISON
za=abs(yy.*exp(j*xx)-z);
k=[];k1=[];i=find(za==min(za));
if ~isempty(i),i1=i(1);  %%NEW
if imag(z(i1))~=0,k=find(z==conj(z(i1)));
if ~isempty(k),k1=k(1);end,end,end %NEW
l=[i1 k1];
if (isempty(k1) & length(i)==1)|length(k)==1 %%NNN
polar(angle(z(l)),abs(z(l)),'oi'),
end %NNN
z(l)=[];end,end

%i=find(za==min(za));k=[];if ~isempty(i),i=i(1);   %%NEW
%if imag(z(i))~=0,k=find(z==conj(z(i)));if ~isempty(k),k=k(1);end,end,end %%NEW
%l=[i k];polar(angle(z(l)),abs(z(l)),'oi'),z(l)=[];end,end

end,end

if bbb==2,title('ADD: LB:POLES  RB:ZEROS  ENTER KEY: PLOTS');
b=0;while b ~= 13
[xx,yy,b]=ginput(1);
if abs(yy)<tol,yy=0;xx=0;end  %%SNAP TO ORIGIN
if abs(yy*sin(xx))<tol,xx=0;end %%SNAP TO REAL VALUE
if xx ~= 0,xx=[xx;-xx];yy=[yy;yy];end
if b==1,polar(xx,yy,'xc5'),pp=yy.*exp(j*xx);p=[p;pp];end
if b==2,if abs(1-abs(yy(1)))<tol,yy=0*yy+1;end,polar(xx,yy,'oc12')
zz=yy.*exp(j*xx);z=[z;zz];end
end,end

if bbb==1 | bbb==2
if isempty(z),num=1;else,num=real(poly(z.'));end                    %%NEW
if isempty(p),den=1;else,den=real(poly(p.'));end                    %%NEW
n=0:100;n0=0*n';n0(1)=1;y=filter(num,den,n0);
w=0:pi/200:pi;zzz=exp(sqrt(-1)*w);h=polyval(num,zzz)./polyval(den,zzz);
hold off,axis('normal');subplot(211)
dtplot(n,y,'.'),grid,title('MENU: LB: DEL  RB: ADD  ENTER KEY: END PROG'),...
xlabel('Impulse response n=0:100'),...
subplot(212),plot(w/2/pi,abs(h)),grid,xlabel('Frequency Response vs F'),...
subplot,hold off,end
[xxx,yyy,bbb]=ginput(1);

end