COMMENT ⊗   VALID 00006 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	begin "AWFIL"
C00003 00003	! requires and other declarations
C00005 00004	! local procedures
C00010 00005	! main program here
C00013 00006	end "AWFIL"
C00014 ENDMK
C⊗;
begin "AWFIL"

require "HEADER.SAI[LIB,ROB]" source_file;

! lotsa stuff drawn from
	JAMLIB.JAM[UP,DOC]
	JAMLIB.SAI[SUB,SYS]
	JAMLIB.HDR[SUB,SYS];

! This program takes the frequency response of a system and calculates the
  coefficients for a set of second order filters which would flatten the
  response of the system.

  The input parameters are:
	a text file named MAGFRQ.TXT which contians a preload_with spec for
	   the magnitude of the frequency array.  Should be NPoints long.

  The outputs are:
	various plots of the system response
	the coefficients of the second order sections

  ;
! requires and other declarations;

require "JAMLIB.REL[SUB,SYS]" library;
    external integer procedure POLY(
	real array A; reference real array ZR,ZI; integer N);
    external procedure RTRAN(
	real array X,A,B; integer log2N; boolean inverse);
    external real procedure AUTSOL(
	real array R,A,K; integer M);
    external string procedure DEFAULT(string Inp; reference string Dev;
		string Devext,Devdev);
    external procedure DSETUP(integer nwds; reference integer id);
    external procedure BUFCLR(integer id,nwds);
    external boolean procedure DRELS(reference integer id);
    external procedure WRITE(integer id,pog);
    external procedure DPYSL(integer id; real array X;
	integer N,Nslices,sliceN; string xlabel,ylabel;
	real vxmin,vxmax);

DEFINE
  NPoints = {512},	comment the number of input points;
  Log2NPoints = {9};

! local procedures;

procedure COMBINE(
    real array Rls,Ims,Ai,Bi; integer M);
    begin "COMBINE"
	boolean array taken[1:M];
	integer ind,upto,jj;

	arrclr(taken);
	ind←1;
	for upto←1 step 1 until M do
	begin "FFCR"
	  if taken[upto] then continue "FFCR";
	  if abs(Ims[upto])<.0000001 then
	  begin "ISR"
	    for jj←upto+1 step 1 until M do
	      if abs(Ims[jj])<.0000001 then done;
	    if jj≤M then
	    begin "ISANO"
		Ai[ind]←-Rls[upto]-Rls[jj];
		Bi[ind]←Rls[upto]*Rls[jj];
		taken[jj]←true;
	    end "ISANO"
	    else begin "NOTAN"
		Ai[ind]←-Rls[upto];
		Bi[ind]←0;
	    end "NOTAN";
	  end "ISR"
	  else begin "ISI"
	    for jj←upto+1 step 1 until M do
	      if abs(Ims[upto]+Ims[jj])<.0000001
		∧ abs(Rls[upto]-Rls[jj])<.0000001 then done;
	    if jj>M then usererr(0,0,"Complex filter?????");
	    Ai[ind]←-2*Rls[upto];
	    Bi[ind]←Rls[upto]↑2+Ims[upto]↑2;
	    taken[jj]←true;
	  end "ISI";
	  ind←ind+1;
	end "FFCR";
    end "COMBINE";

procedure ArrSqr(real array Arr);
  ⊂ "ArrSqr"
  define ndims = {-1}, nwrds = {0}, arlo = {1}, arhi = {2};
  integer LoBound, HiBound, I;
  if arrinfo(Arr,ndims) = 1
    then
      ⊂
      LoBound ← arrinfo(Arr,arlo);
      HiBound ← arrinfo(Arr,arhi);
      for I ← LoBound step 1 until HiBound do Arr[I] ← Arr[I]↑2;
      ⊃
    else
      outstr("ArrSqr: multiply dimensioned array.  No good.  No change."&↓);
  ⊃ "ArrSqr";

procedure GetArray(real array Arr);
  ⊂ "GetArray"
  ! This routine reads in a file (defaults to .FUN extension), which expects
  exactly one "smoothed seg function" as produced by the FUNC program.  The
  function is then loaded into the array Arr;
  integer I,ArrChan,cnt,brk,eof;
  string FilNam,DevNam;
  while true do
    ⊂
    outstr(".FUN file:  ");
    FilNam ← DEFAULT(inchwl, DevNam, "FUN", "DSK");
    outstr("Reading "&DevNam&":"&Filnam&"."&↓);
    open(ArrChan←getchan,DevNam,0,19,19,cnt,brk,eof);
    if eof then errquit(<"Couldn't OPEN file."&↓>);
    lookup(ArrChan,FilNam,eof);
    if ¬eof then done;
    print("File "&DevNam&":"&FilNam&" not found...try again.",↓);
    release(ArrChan);
    ⊃;
    do until realin(ArrChan)= 520.00;
    for I ← 0 step 1 until NPoints do Arr[I] ← realin(ArrChan);
    release(ArrChan);
  ⊃ "GetArray";

procedure Magnitude(real array Ai,Bi,HMag; integer N);
  ⊂ "Magnitude"
  ! Given the coefficients of a cannonical filter (in Ai and Bi), this
    routine calculates the magnitude frequency response, storing the
    result in HMag.  The array HMag is assumed to be N points long.

    |H[i]| = SQRT((.....)/(......))
  ;

  integer OMEGA,I;
  real num1,num2,den1,den2,angle,cospart,sinpart;
  for omega ← 0 step 1 until N do
    begin
    num1←num2←den1←den2←0;
    for I ← 0 step 1 until M do
      begin
      angle ← -I*omega;
      cospart ← cos(agle);
      sinpart ← sin(angle);
      num1 ← num1 + A[I]*cospart;
      num2 ← num2 + A[I]*sinpart;
      den1 ← den1 + B[I]*cospart;
      den2 ← den2 + B[I]*sinpart;
      end;
    H[omega] ← sqrt((num1↑2 + num2↑2)/(den1↑2 + den2↑2));
    end;
  ⊃ "Magnitude";

! main program here;

String str;
integer NCoeffs, brk;
PRINT("How many filter sections:  ");
str ← inchwl;
NCoeffs ← 2*intscan(str,brk);

⊂ "MAIN"

real array MagFrq[0:NPoints];
real array Zeros[0:NPoints];
real array AutCor[1:2*NPoints];
real array APrams[0:NCoeffs],KPrams[1:NCoeffs],TAPrams[0:NCoeffs];
real array RRoots[1:NCoeffs],IRoots[1:NCoeffs];
real array Ai[0:NCoeffs],Bi[0:NCoeffs];
integer id, NRoots, I;

arrclr(Zeros);
DSETUP(2000,id);

GetArray(MagFrq);
DPYSL(id, MagFrq, NPoints, 2, 2, "Frequency", "Magnitude", 0, 512);
WRITE(id,0);
OUTSTR("<alt> for next"&↓);INCHRW;
BUFCLR(id,2000);

! square the magnitude response;
ArrSqr(MagFrq);
DPYSL(id, MagFrq, NPoints, 2, 2, "Frequency", "Magnitude squared", 0, 512);
WRITE(id,0);
OUTSTR("<alt> for next"&↓);INCHRW;
BUFCLR(id,2000);

! take the inverse FFT of the squared magnitude response (yeilds autocorrelation);
RTRAN(AutCor,MagFrq,Zeros,Log2NPoints,TRUE);
DPYSL(id, AutCor, NPoints, 2, 2, "AutCor", "", 0, 512);
WRITE(id,0);
OUTSTR("<alt> for next"&↓);INCHRW;
BUFCLR(id,2000);

! run a Linear Predictor over it to get the filter coefficients;
AUTSOL(AutCor,APrams,KPrams,NCoeffs);

! factor the filter coefficients into their component roots;
for I ← 0 step 1 until NCoeffs do TAPrams[I] ← APrams[NCoeffs-I]; ! invert order;
NRoots ← POLY(TAPrams,RRoots,IRoots,NCoeffs);
If NRoots = 0
  then errquit(<"Couldn't get the roots (?)."&↓>)
  else outstr("Found "&CVS(NRoots)&" roots:"&↓);

! combine the roots into conjugate pairs;
COMBINE(RRoots,IRoots,Ai,Bi,NRoots);

! output the roots, now suitable for a cascaded second order system;
outstr("i"&tab&"Ai"&"Bi"&↓);
for I ← 1 step 1 until NRoots/2 do outstr(I&TAB&cvf(Ai[I])&TAB&cvf(Bi[I])&↓);

⊃ "MAIN";

end "AWFIL";