Hi All,

I bought Fred Harris book Multirate Signal Processing for Communications around 2007 and used it for multirate fir and cic's back then. He had several downloadable programs available for download on the web, one of which is called tony_des.   But I have misplaced them and they don't seem to be available on the web any longer.  Does anyone have those downloads?  In particular I'm interested exploring the recursive all-pass filters in ch10 (using his tony_des2 MATLAB function), but open to additional pointers/suggestions on the topic.  I know that Saramaki has some good papers on this too but I'd like to start with tony_des2.   Much appreciated if someone can provide the downloads for that book.

Thanks,

Chris

Here it is

function tony_des_2(action)
% tony_des_2 forms coefficients of two path recursive polyphase filter
% filter can be halfband quadrature mirror, polynomials in Z^2
% filter can be bandwidth tuned,        Z^-1 => (1-b*z)/(z-b),       b=(1-tan(bw))/(1+tan(bw))
% filter can be center frequency tuned, z^-1 => (-1/z)*(1-cz)/(z-c), c=cos(theta_c)
% filter can be zero-packed             Z^-1 ->  z^-k,               k even     
% order of filter is odd (3,5,7.....)
% wo is normalized band edge frequency .5>wo>.25
%
% %%%%%%%%% this program calls tonycxx.m %%%%%%%%%%%%%
%
% based on paper "Digital Signal Processing Schemes for Efficient Interpolation and Decimation"
% by Valenzuela and Constantinides, IEE Proceedings, Dec 1983
% Script file written by fred harris of SDSU. Copyright 2002

global N_
global ws_
global wb_
global wc_
global kk_
if nargin<1

figure(1)
clf reset
set(1,'name','                Two-Path Allpass Filter Design by fred harris, SDSU','numbertitle','off');

responseplot=axes('units','normalized','position',[.05 .4 .9 .55]);

N_=uicontrol('style','edit','units','normalized','position',[0.18 0.1 0.07 0.07],'string','5');
Ntext=uicontrol('style','text','units','normalized','position',[0.04 0.1 0.13 0.065],'string','Prototype Order (Odd)');

ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string','0.300');
wstext=uicontrol('style','text','units','normalized','position',[0.04 0.17 0.13 0.065],'string','Stopband Edge (0.25-0.5)');

wb_=uicontrol('style','edit','units','normalized','position',[0.42 0.17 0.07 0.07],'string','0.250');
wbtext=uicontrol('style','text','units','normalized','position',[0.28 0.17 0.13 0.065],'string','Passband Edge');

wc_=uicontrol('style','edit','units','normalized','position',[0.66 0.17 0.07 0.07],'string','0.000');
wctext=uicontrol('style','text','units','normalized','position',[0.52 0.17 0.13 0.065],'string','Center Frequency');

kk_=uicontrol('style','edit','units','normalized','position',[0.90 0.17 0.07 0.07],'string','2');
wctext=uicontrol('style','text','units','normalized','position',[0.76 0.17 0.13 0.065],'string','Order (even)');

half_=uicontrol('style','push','units','normalized','position',[0.04 0.25 0.21 0.07],'string','Half-Band Low-Pass','callback','tony_des_2(''start1'')');

tlow_=uicontrol('style','push','units','normalized','position',[0.28 0.25 0.21 0.07],'string','Tune Bandwidth','callback','tony_des_2(''start1'')');

tband_=uicontrol('style','push','units','normalized','position',[0.52 0.25 0.21 0.07],'string','Tune Center Freq','callback','tony_des_2(''start1'')');

fordr_=uicontrol('style','push','units','normalized','position',[0.76 0.25 0.21 0.07],'string','Polynomial Degree','callback','tony_des_2(''start1'')');

tprint=uicontrol('style','push','units','normalized','position',[0.30 0.1 0.65 0.05],'string','Print Coefficients','callback','tony_des_2(''start2'')');

readtext=uicontrol('style','text','units','normalized','position',[0.280 0.03 0.65 0.05],'string','(Highlight Parameter Value, Enter Desired Value, Depress Button)');
action='start1';
end

if strcmp(action,'start1')
   coef=0;
end
if strcmp(action,'start2')
   coef=1;
end

      N=str2num(get(N_,'string'));
   if rem(N,2)~=1
      N=N+1;
      N_=uicontrol('style','edit','units','normalized','position',[0.18 0.1 0.07 0.07],'string',N);
   end
   ws=str2num(get(ws_,'string'));
   if ws<0.25
      ws=0.300;
      ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string',ws);
   end
   if ws>0.5
      ws=0.300;
      ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string',ws);
   end
   wb=str2num(get(wb_,'string'));
   if wb>0.5
      wb=0.25;
      wb_=uicontrol('style','edit','units','normalized','position',[0.42 0.17 0.07 0.07],'string','0.250');
   end
   wc=str2num(get(wc_,'string'));
   if wc>0.5
      wc=0.0;
      wc_=uicontrol('style','edit','units','normalized','position',[0.66 0.17 0.07 0.07],'string','0.000');
   end
  
   kk=str2num(get(kk_,'string'));
   if rem(kk,2)~=0;
      kk=2;
      kk_=uicontrol('style','edit','units','normalized','position',[0.90 0.17 0.07 0.07],'string','2');
   end

   if wb==0.25 & wc==0
      mode=1;
   else
      mode=2;
   end
  
   aa=1;
   if wc==0;
      aa=-aa;
   end
  
[rts1,rts2,coef0,coef1,ff1,ff2]=tonycxx(N,ws,wb,wc,kk,mode);

figure(2)
plot(roots(rts1),'x','linewidth',2,'markersize',8);
hold on;
plot(roots(rts2),'x','linewidth',2,'markersize',8);
plot(roots(conv(rts1,fliplr(rts2))-aa*conv(fliplr(rts1),rts2)),'o','linewidth',2)
plot(exp(j*2*pi*[0:.01:1]),'r');
plot([-1.1 1.1],[0 0],'k')
plot([0 0],[-1.1 1.1],'k')

hold off;
grid
axis([-1.2 1.2 -1.2 1.2]);
axis('square')
title('Roots of Two-Path Filter (Sum of Paths)')

figure(1)
plot(0:1/1024:.5-1/1024,20*log10(abs(ff1+0.0000001)));
hold
plot(0:1/1024:.5-1/1024,20*log10(abs(ff2+0.0000001)),'r--');
axis([0 .5 -100 10]);

hold

grid;
title('Magnitude Response of Two Path Filter');
xlabel('Normalized Frequency (f/fs)');
ylabel('Log Magnitude (dB)');

% write out coefficients
format long
if coef==1
    disp(' ')
    disp('Denominator Polynomials, path-0')
        if wc==0
        disp('   (a0 Z^2  +  a1 Z^1  +  a2)')
        else
        disp('   (a0 Z^4  +  a1 Z^3  +  a2 Z^2  +  a3 Z^1  +  a4)')
end
   disp(' ')
   disp(coef0)
   disp(' ')
  
    disp('Denominator Polynomials, path-1')
    if wc==0
    disp('    (b0 Z^2  +  b1 Z^1  +  b2)')
    else
    disp('    (b0 Z^4  +  b1 Z^3  +  b2 Z^2  +  b3 Z^1  +  b4)')
end
   disp(' ')  
   disp(coef1)
end

Thanks!  Looks like it also needs the tonycxx() function.

function [rts1,rts2,coef0,coef1,ff1,ff2]=tonycxx(n_ord,wo,bw,fc,k_ordr, mode)
% mode = 1 for lowpass prototype
% mode = 2 to bandwidth and center frequency of prototype
% k_ordr for filter order
% function called by tony_des_2,
% computes paramters of two-path recursive all-pass filter
% also computes parameters when frequency transformed to two-path
% low-pass or band-pass filter
% Script file written by fred harris of SDSU. Copyright 2002

mode;
wt=4*pi*(wo-0.25);
k=tan(0.25*(pi-wt));
k=k*k;

kk=sqrt(1-k*k);

e=0.5*(1-sqrt(kk))/(1+sqrt(kk));

q=e+2*(e^5)+15*(e^9)+150*(e^13);
ww=zeros(1,(n_ord-1)/2);
aa=ww;
cc=[ww 0];
%step 2

for ii=1:(n_ord-1)/2

ww(ii)=2*(q^0.25)*(sin(pi*ii/n_ord)-(q^2)*sin(3*pi*ii/n_ord));
ww(ii)=ww(ii)/(1-2*(q*cos(2*pi*ii/n_ord)-(q^4)*cos(4*pi*ii/n_ord)));

wwsq=ww(ii)*ww(ii);
aa(ii)=sqrt(((1-wwsq*k)*(1-wwsq/k)))/(1+wwsq);
cc(ii)=(1-aa(ii))/(1+aa(ii));
end

ordr1=floor((n_ord-1)/4);
ordr2=ordr1;
if n_ord-1-4*ordr1 ~= 0
 ordr1=ordr1+1;
end

coef0=zeros(ordr1,3);
coef1=zeros(ordr2,3);

zz=zeros(1,k_ordr-1);

for ii=1:ordr1
den0(ii,:)=[1 zz cc(2*ii-1)];
end
coef0=den0;

for ii=1:ordr2
den1(ii,:)=[1 zz cc(2*ii)];
end
coef1=den1;

h0=[1];
for ii=1:ordr1
h0=conv(h0,den0(ii,:));
end

h1=[1];
for ii=1:ordr2
h1=conv(h1,den1(ii,:));
end

zz2=zeros(1,k_ordr/2);
h1=[h1 zz2];

g0=fliplr(h0);
g1=fliplr(h1);

rts1=h0;
rts2=h1;

[tp,ww]=freqz(g0,h0,512);
[bt,ww]=freqz(g1,h1,512);

ff1=0.5*(tp+bt);
ff2=0.5*(tp-bt);

if mode==1
return  

else

tt=tan(pi*bw);
b=(1-tt)/(1+tt);

c=cos(2*pi*fc);

if fc==0
  
den0=zeros(ordr1,k_ordr+1);
den1=zeros(ordr2,k_ordr+1);
zz=zeros(1,(k_ordr/2)-1);
   for nn=1:ordr1
c00=1+cc(2*nn-1)*b*b;
c01=-2*b*(1+cc(2*nn-1));
c01=c01/c00;
c02=cc(2*nn-1)+b*b;
c02=c02/c00;
den0(nn,:)=[1 zz c01 zz c02];
   end
coef0=den0;
  
   for nn=1:ordr2
c10=1+cc(2*nn)*b*b;
c11=-2*b*(1+cc(2*nn));
c11=c11/c10;
c12=cc(2*nn)+b*b;
c12=c12/c10;
den1(nn,:)=[1 zz c11 zz c12];
end
zz2=zeros(1,k_ordr/2-1);
den1(ordr2+1,1:1+k_ordr/2)=[1 zz2 -b];  
coef1=den1;

h0=[1];
for ii=1:ordr1
   h0=conv(h0,den0(ii,:));
end

h1=[1 zz2 -b];
for ii=1:ordr2
   h1=conv(h1,den1(ii,:));
end

g0=fliplr(h0);
g1=fliplr(h1);
rts1=h0;
rts2=h1;

[tp,ww]=freqz(g0,h0,512);
[bt,ww]=freqz(g1,h1,512);

ff1=0.5*(tp+bt);
ff2=0.5*(tp-bt);

else
  
den0=zeros(ordr1,2*k_ordr+1);
den1=zeros(ordr2,2*k_ordr+1);

zz=zeros(1,(k_ordr/2)-1);

for nn=1:ordr1
c00=1+cc(2*nn-1)*b*b;
c01=-2*c*(1+b)*(1+cc(2*nn-1)*b);
c01=c01/c00;
c02=(1+cc(2*nn-1))*(c*c*(1+b*b)+2*b*(1+c*c));
c02=c02/c00;
c03=-2*c*(1+b)*(cc(2*nn-1)+b);
c03=c03/c00;
c04=cc(2*nn-1)+b*b;
c04=c04/c00;
den0(nn,:)=[1 zz c01 zz c02 zz c03 zz c04];
end
coef0=den0;

for nn=1:ordr2
c10=1+cc(2*nn)*b*b;
c11=-2*c*(1+b)*(1+cc(2*nn)*b);
c11=c11/c10;
c12=(1+cc(2*nn))*(c*c*(1+b*b)+2*b*(1+c*c));
c12=c12/c10;
c13=-2*c*(1+b)*(cc(2*nn)+b);
c13=c13/c10;
c14=cc(2*nn)+b*b;
c14=c14/c10;
den1(nn,:)=[1 zz c11 zz c12 zz c13 zz c14];
end
zz2=zeros(1,k_ordr/2-1);
den1(ordr2+1,1:1+k_ordr)=[1 zz2 -c*(1+b) zz2 b];
coef1=den1;

h0=[1];
for ii=1:ordr1
   h0=conv(h0,den0(ii,:));
end

h1=[1 zz -c*(1+b) zz b];
for ii=1:ordr2
   h1=conv(h1,den1(ii,:));
end

g0=fliplr(h0);
g1=fliplr(h1);
rts1=h0;
rts2=h1;

[tp,ww]=freqz(g0,h0,512);
[bt,ww]=freqz(g1,h1,512);

ff1=0.5*(tp-bt);
ff2=0.5*(tp+bt);
end

end

Got it - Thanks!