From 22f30687ec3544cd7d27d49b05b2bf8953574385 Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Wed, 10 Jan 2024 09:43:53 +0100
Subject: [PATCH] All pass IIR filter desing by F. Harris.
---
applications/lofar2/model/harris/tony_des_2.m | 160 ++++++++++++++
applications/lofar2/model/harris/tonycxx.m | 201 ++++++++++++++++++
2 files changed, 361 insertions(+)
create mode 100644 applications/lofar2/model/harris/tony_des_2.m
create mode 100644 applications/lofar2/model/harris/tonycxx.m
diff --git a/applications/lofar2/model/harris/tony_des_2.m b/applications/lofar2/model/harris/tony_des_2.m
new file mode 100644
index 0000000000..3789844d63
--- /dev/null
+++ b/applications/lofar2/model/harris/tony_des_2.m
@@ -0,0 +1,160 @@
+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
diff --git a/applications/lofar2/model/harris/tonycxx.m b/applications/lofar2/model/harris/tonycxx.m
new file mode 100644
index 0000000000..3d89ffd3e3
--- /dev/null
+++ b/applications/lofar2/model/harris/tonycxx.m
@@ -0,0 +1,201 @@
+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
--
GitLab