All pass IIR filter desing by F. Harris.

applications/lofar2/model/harris/tony_des_2.m

+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

applications/lofar2/model/harris/tonycxx.m

+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