gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\maxflat.m
function [b,a,b1,b2] = maxflat(Z,P,wo,PF)
%MAXFLAT Maximally flat (a.k.a. generalized Butterworth) digital filter design.
% [B,A] = MAXFLAT(NB,NA,Wn) is a lowpass Butterworth filter with numerator
% and denominator coefficients B and A of orders NB and NA respectively.
% Wn is the cut-off frequency at which the filter's magnitude response
% is equal to 1/sqrt(2) (approx. -3 dB). Wn must be between 0 and 1, with
% 1 corresponding to half the sample rate of the filter. When NA=0, the
% range of Wn is further restricted due to the smoothness of the filter's
% response.
%
% B = MAXFLAT(NB,'sym',Wn) is a symmetric FIR Butterworth filter. NB must be
% even, and Wn is restricted to a subinterval of [0,1]. The function will
% raise an error if you specify Wn outside of this subinterval.
%
% [B,A,B1,B2] = MAXFLAT(NB,NA,Wn) and [B,A,B1,B2] = MAXFLAT(NB,'sym',Wn)
% return two polynomials B1 and B2 whose product is equal to the numerator
% polynomial B (i.e., B = CONV(B1,B2)). B1 contains all the zeros at z=-1,
% and B2 contains all the other zeros.
%
% EXAMPLES
% NB = 10; NA = 2; Wn = 0.6;
% [b,a,b1,b2] = maxflat(NB,NA,Wn);
% h = maxflat(NB,'sym',Wn);
%
% MONITORING THE DESIGN
% For a textual display of the design table used in the design, use a
% trailing 'trace' argument, e.g. MAXFLAT(NB,NA,Wn,'trace'). To display
% plots of the filter's magnitude, group delay and zeros and poles, use
% a trailing 'plots' argument, MAXFLAT(NB,NA,Wn,'plots'). For both text
% and plots, use 'both'.
%
% See also BUTTER, FREQZ, FILTER.
% by I. W. Selesnick and C. S. Burrus
% see: Generalized Digital Butterworth Filter Design
% by I. W. Selesnick and C. S. Burrus,
% IEEE International Conference on Acoustics,
% Speech and Signal Processing,
% May 1996, Volume 3, p. 1367.
% Author: Ivan Selesnick, Rice University, September 1995.
% Copyright (c) 1988-98 by The MathWorks, Inc.
% $Revision: 1.1 $ $Date: 1998/06/03 14:43:16 $
% required subprograms : spec_table.m, choose.m
TEXT_PF = 0;
PLOT_PF = 0;
SYMM = 0;
if isstr(P)
SYMM = 1;
P = 0;
if rem(Z,2)==1
error('For symmetric FIR filters, only even orders are allowed.')
else % (Z is even)
Z = Z/2;
end
end
if nargin > 3
if ~isstr(PF)
error('Plot flag must be ''trace'', ''plots'', or ''both''.')
end
switch lower(PF(1))
case 't'
TEXT_PF = 1;
case 'p'
PLOT_PF = 1;
case 'b'
TEXT_PF = 1;
PLOT_PF = 1;
otherwise
error('Plot flag must be ''trace'', ''plots'', or ''both''.')
end
end
SN = 1e-7;
table = spec_table(Z,P,SYMM);
wo = wo*pi;
k = max(find((table(:,4) < wo/pi+SN) & (table(:,5) > wo/pi-SN)));
%%%%%%%%%%%%%%%%%%%% DISPLAY TABLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if TEXT_PF
disp(' Table:')
disp(' ')
disp(' L M N wo_min/pi wo_max/pi')
disp(' ')
disp(table)
if SYMM
disp(' L and M are doubled for symmetric FIR filters.')
disp(' ')
end
% L+M : total number of zeros
% N : total number of poles
%
% L : number of zeros at z=-1
% M : number of zeros contributing to passband flatness
% N : number of poles
end
if PLOT_PF
clf
set(gcf,'nextplot','add');
a1 = axes('units','normalized','position',[.1 .6 .8 .3]);
a2 = axes('units','normalized','position',[.1 .1 .35 .4]);
a3 = axes('units','normalized','position',[.55 .1 .35 .4]);
end
%%%%%%%%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(k) == 0
error(sprintf('For NB=%g and NA=%g, Wn must be between %g and %g.',...
(1+SYMM)*Z,P,table(1,4),table(end,5)))
end
%%%%%%%%%%%%%%%%%%%% CALCULATE x-Domain function %%%%%%%%%%%%%%%%%%
L = table(k,1);
M = table(k,2);
N = table(k,3);
xo = (1-cos(wo))/2;
%Fo = (1/2)^2;
Fo = 1/2;
if SYMM
Fo = sqrt(Fo);
end
if N == 0
%%%%%%%%%%%% FIR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = M-1:-1:0;
R = [0, choose(M-k-1,0).*choose(L+k-1,k)];
T = [choose(M-k-2,-1).*choose(L+k,k), 0];
c1 = (Fo/(1-xo)^L - polyval(R,xo));
c2 = polyval(T,xo);
S = c2*R + c1*T;
Q = 1;
elseif M == 0
%%%%%%%%%%%% ALL ZEROS at z = -1 %%%%%%%%%%%%%%%%%%%%%%%%%%
k = min([L,N]):-1:0;
Q = choose(L,k).*(-1).^k;
if L < N
Q = [zeros(1,N-L), Q];
end
c1 = ((1/Fo)*(1-xo)^L - polyval(Q,xo))/(xo^N);
Q(1) = Q(1) + c1;
S = 1;
else % M > 0
%%%%%%%%%%%% Some Zeros contribute to passband %%%%%%%%%%%%
k = M-1:-1:0;
R = [0, choose(M+N-k-1,N).*choose(L-N+k-1,k)];
T = [choose(M+N-k-2,N-1).*choose(L-N+k,k), 0];
k = N:-1:0;
AGR = choose(M+N-k-1,M-1).*choose(M+L-1,k).*(-1).^k;
AGT = choose(M+N-k-1,M-1).*choose(M+L-1,k-1).*(-1).^(k-1);
c1 = (((1-xo)^L)*polyval(R,xo) - Fo*polyval(AGR,xo));
c2 = (Fo*polyval(AGT,xo) - ((1-xo)^L)*polyval(T,xo));
S = c2*R + c1*T;
Q = c2*AGR + c1*AGT;
end
%%%%%%%%%%%%%%%%%%%% convert to z-domain %%%%%%%%%%%%%%%%%%%%%%%%%%
tmp = 1 - 2*roots(S);
br1 = tmp+sqrt(tmp.^2-1);
br2 = tmp-sqrt(tmp.^2-1);
br = sort([br1; br2]+eps*sqrt(-1)); % sort by absolute value
if ~SYMM
if c1 ~= 0
br = br(1:M); % take roots INSIDE unit circle
else
br = [br(1:M-1); 0]; % take roots INSIDE unit circle
end
end
b2 = real(poly(br));
Qrts = roots(Q);
ar1 = 1-2*Qrts+sqrt((1-2*Qrts).^2-1);
ar2 = 1-2*Qrts-sqrt((1-2*Qrts).^2-1);
ar = sort([ar1; ar2]+eps*sqrt(-1)); % sort by absolute value
if ~SYMM
ar = ar(1:N); % take roots INSIDE unit circle
end
a = real(poly(ar));
%%%%%%%%%%%%%%%%%%%% construct b1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SYMM
L = 2*L;
end
b1 = 1;
b = b2;
for t = 1:L
b1 = conv(b1,[1 1]);
b = conv(b,[1 1]);
end
%%%%%%%%%%%%%%%%%%%% normalize %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if abs(sum(b2)) < 100*eps
% terminate program
error('A pole-zero pair cancels at z=1. Use a lower order.')
else
b2 = b2*sum(a)/(sum(b1)*sum(b2));
b = conv(b1,b2);
end
%%%%%%%%%%%%%%%%%%%% DISPLAY RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(a) & PLOT_PF
% ---- PLOT MAGNITUDE ----
[Mb1,w] = h2mag(b1);
[Mb2,w] = h2mag(b2);
[Ma,w] = h2mag(a);
Mag = Mb1.*Mb2./Ma;
axes(a1)
plot(w/pi,Mag)
xlabel('\omega/\pi')
ylabel('Magnitude')
title('Frequency response')
axis([0 1 -.1 1.1])
% ---- PLOT POLE-ZERO DIAGRAM ----
I = sqrt(-1);
SN = 1e-7;
axes(a2)
circ = exp([0:100]/100 * 2 * pi * I);
plot(circ,'--')
hold on
if length(ar)
plot(ar+SN*I,'x')
end
if length(br)
plot(br+SN*I,'o')
end
plot([-1+SN*I],'o')
hold off
axis([-1 1 -1 1]*1.2)
xlabel('Real')
ylabel('Imaginary')
if SYMM
title('Zero plot (zeros outside circle not shown)')
else
title('Pole-zero plot')
end
s1 = ['<- deg ',num2str(L)];
t1 = text(-0.9,0,s1);
axis('square')
% ---- PLOT GROUP DELAY ----
gd = grpdelay(b2,a,w)+L/2;
axes(a3)
plot(w/pi,gd)
mingd = min(gd);
maxgd = max(gd);
if (maxgd-mingd)<(1e-10)*(abs(maxgd)+abs(mingd))
set(a3,'ylim',[0 2*maxgd])
end
xlabel('\omega/\pi')
ylabel('Samples')
title('Group delay')
set(gcf,'nextplot','replace');
end
function table = spec_table(Z,P,SYMM)
%
% table = spec_table(Z,P,SYMM);
%
% General Butterworth Filter Design
% This program decides how many zeros should lie at
% z = -1 and how many should contribute to the passband
%
% Z : total number of zeros
% P : total number of poles
%
% % Example:
%
% Z = 10; P = 2;
% table = spec_table(Z,P);
%
% Ivan Selesnick, Rice University,
% September 1995.
% see: Generalized Digital Butterworth Filter Design
% required subprograms : choose.m
if Z <= P
table = [Z 0 P 0 1];
return
else
table = [];
end
%Fo = (1/2)^2; % value of mag squared at half-mag freq.
Fo = 1/2; % value of mag squared at half-mag freq.
if SYMM
Fo = sqrt(Fo);
end
%%%%%%% begin with all zeros at z = -1 %%%%%%%
L = Z;
M = 0;
N = P;
if N > 0
if rem(N,2) == 0
c = 0; % N even
else
c = choose(L-1,N); % N odd
end
%%%%%%% construct polynomial for checking boundary %%%%%%
Y = zeros(1,L+1);
for i = 0:N
Y(i+1) = choose(L,i)*(1-Fo)*(-1)^i;
end
for i = N+1:L
Y(i+1) = choose(L,i)*(-1)^i;
end
Y(N+1) = Y(N+1) - c*Fo;
Y = Y(L+1:-1:1);
%%% extract appropriate root %%%
r = specroot(Y);
wo_max = acos(1-2*r);
table = [Z, 0, P, 0, wo_max/pi];
end
%%%%%%%%% go through each and check %%%%%%%%%
for M = 1:Z-P-1
L = Z-M;
k = M-1:-1:0;
R = [0, choose(M+N-k-1,N).*choose(L-N+k-1,k)];
T = [choose(M+N-k-2,N-1).*choose(L-N+k,k), 0];
k = M+L:-1:0;
GR = choose(M+N-k-1,M-1).*choose(M+L-1,k).*(-1).^k;
GT = choose(M+N-k-1,M-1).*choose(M+L-1,k-1).*(-1).^(k-1);
AGR = GR(M+L+1-N:M+L+1);
AGT = GT(M+L+1-N:M+L+1);
%%% ------- wo_max ---------------------------------- %%%
if rem(N,2) == 0
c = (L-N)/(M+N); % N even
else
c = (L-N)/N; % N odd
end
Y = [zeros(1,M+L-N), Fo*AGR+c*Fo*AGT]-GR-c*GT;
%%%%%%% extract appropriate root %%%
r = specroot(Y);
wo_max = acos(1-2*r);
%%% ------- wo_min ---------------------------------- %%%
if rem(N,2) == 0
c = -1; % N even
Y = [zeros(1,M+L-N), Fo*AGR+c*Fo*AGT]-GR-c*GT;
else
c = inf; % N odd
Y = [zeros(1,M+L-N), Fo*AGT]-GT;
end
%%%%%%% extract appropriate root %%%
r = specroot(Y);
wo_min = acos(1-2*r);
table = [table; [L,M,N, wo_min/pi, wo_max/pi]];
end
if N > 0
M = Z-P;
L = Z-M;
table = [table; [L,M,N, table(Z-P,5), 1]];
end
function a = choose(n,k)
%
% a = choose(n,k)
% BINOMIAL COEFFICIENTS
%
% allowable inputs:
% n : integer, k : integer
% n : integer vector, k : integer
% n : integer, k : integer vector
% n : integer vector, k : integer vector (of equal dimension)
%
nv = n;
kv = k;
if (length(nv) == 1) & (length(kv) > 1)
nv = nv * ones(size(kv));
elseif (length(nv) > 1) & (length(kv) == 1)
kv = kv * ones(size(nv));
end
a = nv;
for i = 1:length(nv)
n = nv(i);
k = kv(i);
if n >= 0
if k >= 0
if n >= k
c = prod(1:n)/(prod(1:k)*prod(1:n-k));
else
c = 0;
end
else
c = 0;
end
else
if k >= 0
c = (-1)^k * prod(1:k-n-1)/(prod(1:k)*prod(1:-n-1));
else
if n >= k
c = (-1)^(n-k)*prod(1:-k-1)/(prod(1:n-k)*prod(1:-n-1));
else
c = 0;
end
end
end
a(i) = c;
end
function r = specroot(Y)
r = roots(Y);
r = r(imag(r)==0);
[temp,i] = min(abs(r-0.5));
r = r(i);
% other alternatives:
%SN = 1e-5;
%r = newton(Y,[SN 1-SN]);
%f = inline('polyval(p,x)','x','p');
%r = fzero(f,.5,eps,0,Y);
function [M,w] = h2mag(h)
% [M,w] = h2mag(h)
L = 2^9;
f = fft(h,2*L);
M = abs(f);
M = M(1:L+1);
w = [0:L]'*pi/L;