apeqp2.m

分数延时的matlab代码

function ap = apeqp2(Nap,x,wp)
% APEQP2.M
% MATLAB m-file for ALLPASS fractional delay design
% using EQUIRIPPLE PHASE approximation
% format: ap = apeqp2(Nap,x,wp)
% Input: Nap = filter order N
%        x = fractional delay (0 < x <= 1) 
%        wp = passband edge of approximation (0 < wp < 1)      
% Output: Filter coefficient vector ap = ap(1)...ap(N+1)
%         (Note that ap(1)=1 always)
% Subroutines: standard MATLAB & SP Toolbox functions &
%              numint2.m, eigsolv2.m, envelop2.m
% Timo Laakso 27.02.1993
% Revised 15.01.1996 by Timo Laakso
%         17.01.1996 by Vesa Valimaki

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Initialization: 
%
dw = wp/200;         % 201-point approximation grid
%dw = wp/511;         % 512-point approximation grid (VPV 9.12.1996)
wf = [0:dw:wp];      % vector of freq. pts
%
thP = -(Nap+x)*wf;  % prescribed phase
%
% beta and initializations
%
beta = -0.5*x*wf; beta2=2*beta;
sinB = sin(beta); cosB = cos(beta);
wc = wf-wf; WW = wc; W0 = wc+1-wc;     % weight vectors
ap = zeros(1,Nap+1); ap(1)=1;          % allpass coefficient vector
P4 = zeros(Nap+1);                     % eigenfilter matrix
Nap1=Nap-1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute P4:
%
EP = 2;     % initial exponent parameter     
% initial value ap from P1!
lser=1000; it=0;
%while ((lser>1.0E-8)&(it<51));
while ((lser>1.0E-12)&(it<171));  % VPV 9.12.1996
  if (EP<51)
    if (lser<1.0E-4)
      if (EP<4) EP = EP+1;
      elseif (EP<10) EP = EP+2;
      else EP = EP+5;
      end;      % 3rd if
      lser = 1000;
    end         % 2nd if
  end           % 1st if
  it = it+1;  % iteration number 
  Hap = freqz(ap,1,wf);   % weight function from the previous solution
  wc = cosB.*real(Hap) - sinB.*imag(Hap);
  thAP4 = -Nap*wf+2*unwrap(angle(freqz(1,ap,wf)));
  deltaT4 = thP-thAP4;
  wenv = envelop2(deltaT4); % create the envelope for equiripple weighting
  W = (wenv.^EP)./(wc.^2);
  for k=1:Nap+1
    k1wf = (k-1)*wf;
    cT(k) = numint2(W.*cos(k1wf)); 
    cH(k) = numint2(W.*cos(beta2-k1wf));
    rH(k) = numint2(W.*cos(beta2-(k+Nap1)*wf));
  end % for k
  T1 = toeplitz(cT); H1 = hankel(cH,rH);
  P4 = T1-H1; apold = ap;
  ap= eigsolv2(P4,Nap+1);     % solve the eigenvalue problem
  ap = (ap+(EP-1)*apold)/EP;
  dif=apold-ap; lser=dif*dif';
end              % while lser,it
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%