<*+M2EXTENSIONS *>
<*-CHECKDIV *>
<*-CHECKRANGE *>
<*-COVERFLOW *>
<*-IOVERFLOW*>

<*+NOPTRALIAS*>
<*CPU="PENTIUMPRO"*>
<*-DOREORDER *>
<*-CHECKNIL *>
<*-CHECKSET*>
<*-PROCINLINE*>
<* IF __GEN_C__ THEN *>
  <*-GENCTYPES*>
  <*+COMMENT*>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-GENDATE*>
  <*-LINENO*>
  <*-CHECKINDEX*>
  <*+GENCDIV*>
  <*-GENKRC*>
  <*+NOOPTIMIZE*>
  <*-GENSIZE*>
  <*-ASSERT*>
<* ELSE *>
  <*-GENHISTORY*>
  <*-GENDEBUG*>
  <*-LINENO*>
  <*-CHECKINDEX*>
<* END *>

IMPLEMENTATION MODULE sdr; (* rtl_tcp iq fm demodulator by OE5DXL *)


FROM tcpb IMPORT connecttob, readsockb;
FROM tcp IMPORT sendsock;
FROM SYSTEM IMPORT INT8, CARD8, INT16, CARD16, CAST, ADR;
FROM osi IMPORT OpenWrite, OpenRead, WrBin, Close, sqrt, ln, pi, sin, cos,
     WrFixed, WrStrLn, WrInt, WrStr, usleep, Werr, ALLOCATE, DEALLOCATE;
FROM aprsstr IMPORT Assign;
FROM unistd IMPORT read;

CONST
      IQBUF=65536;
      DDSMAXLEN=2048;
      AFCSPEED=400000000;                                 (* slower afc reaction *)
      AFCRECENTER=1024;                                   (* pull afc to middle *)
      SSBDDSSIZE=2048;
      SSBDDSSIZE4=SSBDDSSIZE DIV 4;
      AFCRSSICLAMP=0.7;                                   (* until rssi level below peak to do afc *)
      AFCRSSISPEED=0.9985;                                (* afc peak rssi fall speed *)
      PI=3.1415926535;
      PI2=PI*2.0;
      IQLEVEL=1.0/64.0;                                   (* correct iq level on input format *)

TYPE 
     SET32=SET OF [0..31];
     Complex=RECORD
               Re,
               Im: REAL;
             END;
     Complexu8=RECORD
               Re,
               Im:CARD8;
             END;
     Complexi16=RECORD
               Re,
               Im: INT16;
             END;
     COMPI=RECORD
               r,
               i: INTEGER;
             END;


VAR 
    fd:INTEGER;
    isfile:BOOLEAN;
    audiohz,
    rtlhz, bytespersamp: CARDINAL;
    reconnect:BOOLEAN;
    iqbuf  :RECORD CASE :CARDINAL OF
                       0:b:ARRAY[0..IQBUF-1] OF CHAR;
                      |1:u8:ARRAY[0..IQBUF DIV 2-1] OF Complexu8;
                      |2:i16:ARRAY[0..IQBUF DIV 2-1] OF Complexi16;
                      |3:i32:ARRAY[0..IQBUF DIV 2-1] OF COMPI;
                      |4:f32:ARRAY[0..IQBUF DIV 2-1] OF Complex;
                     END;
              END;
    sndw:CARDINAL;
    DDS: ARRAY[0..DDSMAXLEN-1] OF INT16;
    DDSR: ARRAY[0..DDSMAXLEN-1] OF REAL;
    url:ARRAY[0..1000] OF CHAR;
    port:ARRAY[0..10] OF CHAR;
    ddslen:SET32;
    ddslen4:CARDINAL;
    SSBDDS:ARRAY[0..SSBDDSSIZE-1] OF REAL;
    onesamp,
    sampfrac,
    iqlevel,
    varbwk:REAL;

PROCEDURE initdds(size:CARDINAL);
VAR i:CARDINAL;
    d:REAL;
BEGIN
  IF size>HIGH(DDS) THEN size:=HIGH(DDS)+1 END;
  d:=2.0*pi/FLOAT(size);
  FOR i:=0 TO size-1 DO
    DDSR[i]:=32767.5*sin(FLOAT(i)*d);
    DDS[i]:=VAL(INTEGER, DDSR[i]);
  END;
  ddslen:=CAST(SET32, size-1);
  ddslen4:=size DIV 4;
END initdds;

PROCEDURE initssbdds(VAR dds:ARRAY OF REAL);
VAR i:CARDINAL;
    d:REAL;
BEGIN
  d:=2.0*pi/FLOAT(HIGH(dds)+1);
  FOR i:=0 TO HIGH(dds) DO dds[i]:=sin(FLOAT(i)*d) END;
END initssbdds;


PROCEDURE iir512(rx:pRX; a,b:CARDINAL);
CONST FG=-9;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;

BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);

  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir512;


PROCEDURE iir256(rx:pRX; a,b:CARDINAL);
CONST FG=-8;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);
  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir256;


PROCEDURE iir128(rx:pRX; a,b:CARDINAL);
CONST FG=-7;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);

  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir128;

PROCEDURE iir64(rx:pRX; a,b:CARDINAL);
CONST FG=-6;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);

  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir64;

PROCEDURE iir32(rx:pRX; a,b:CARDINAL);
CONST FG=-5;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);

  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir32;


PROCEDURE iir16(rx:pRX; a,b:CARDINAL);
CONST FG=-4;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);
  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir16;


PROCEDURE iir8(rx:pRX; a,b:CARDINAL);
CONST FG=-3;
VAR i,ph,dfc:CARDINAL;
    ddsr, ddsi, r1, r2, r3, i1, i2, i3:INTEGER;
    x:COMPI;
BEGIN
  WITH rx^.tapre DO r1:=uc1; r2:=uc2; r3:=il END;
  WITH rx^.tapim DO i1:=uc1; i2:=uc2; i3:=il END;
  ph:=rx^.phase;
  dfc:=CAST(CARDINAL, CAST(SET32, rx^.df+CAST(CARDINAL, rx^.afckhz))*ddslen);
  i:=a;
  REPEAT
    x:=iqbuf.i32[i];
    INC(i);

    ddsr:=DDS[ph];                                               (* mix osz sin *)
    ddsi:=DDS[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];   (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + ASH(x.r*ddsr - x.i*ddsi -r1-r3, FG);                  (* mixer + lowpass i *)
    r2:=r2 + ASH(r3-r2, FG);
    r3:=r3 + ASH(r1-r2, FG+1);
    i1:=i1 + ASH(x.r*ddsi + x.i*ddsr -i1-i3, FG);                  (* mixer + lowpass q *)
    i2:=i2 + ASH(i3-i2, FG);
    i3:=i3 + ASH(i1-i2, FG+1);

  UNTIL i>=b;
  WITH rx^.tapre DO uc1:=r1; uc2:=r2; il:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO uc1:=i1; uc2:=i2; il:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iir8;


PROCEDURE iirvar(rx:pRX; a,b:CARDINAL; bw:REAL);                 (* variable bandwidth *)
CONST FG=128;
VAR i,ph,dfc:CARDINAL;
    xr, xi, r1, r2, r3, i1, i2, i3, bw2, ddsr, ddsi:REAL;
BEGIN
  WITH rx^.tapre DO r1:=ucr1; r2:=ucr2; r3:=ilr END;
  WITH rx^.tapim DO i1:=ucr1; i2:=ucr2; i3:=ilr END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  bw2:=bw*2.0;
  REPEAT
    xr:=VAL(REAL, iqbuf.i32[i].r);
    xi:=VAL(REAL, iqbuf.i32[i].i);
    INC(i);

    ddsr:=DDSR[ph];                                              (* mix osz sin *)
    ddsi:=DDSR[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];  (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + (xr*ddsr - xi*ddsi -r1-r3)*bw;                      (* mixer + lowpass i *)
    r2:=r2 + (r3-r2)*bw;
    r3:=r3 + (r1-r2)*bw2;
    i1:=i1 + (xr*ddsi + xi*ddsr -i1-i3)*bw;                      (* mixer + lowpass q *)
    i2:=i2 + (i3-i2)*bw;
    i3:=i3 + (i1-i2)*bw2;

  UNTIL i>=b;
  WITH rx^.tapre DO ucr1:=r1; ucr2:=r2; ilr:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO ucr1:=i1; ucr2:=i2; ilr:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iirvar;

PROCEDURE iirvarf4(rx:pRX; a,b:CARDINAL; bw:REAL);                 (* variable bandwidth *)
CONST FG=128;
VAR i,ph,dfc:CARDINAL;
    r1, r2, r3, i1, i2, i3, bw2, ddsr, ddsi:REAL;
    x:Complex; 
BEGIN
  WITH rx^.tapre DO r1:=ucr1; r2:=ucr2; r3:=ilr END;
  WITH rx^.tapim DO i1:=ucr1; i2:=ucr2; i3:=ilr END;
  ph:=rx^.phase;
  dfc:=rx^.df+VAL(CARDINAL, rx^.afckhz);
  i:=a;
  bw2:=bw*2.0;
  REPEAT
    x:=iqbuf.f32[i];
    INC(i);

    ddsr:=DDSR[ph];                                              (* mix osz sin *)
    ddsi:=DDSR[CAST(CARDINAL, CAST(SET32, ph+ddslen4)*ddslen)];  (* mix osz cos *)
    ph:=CAST(CARDINAL, CAST(SET32, ph+dfc)*ddslen);              (* drive mix osz *)

    r1:=r1 + (x.Re*ddsr - x.Im*ddsi -r1-r3)*bw;                      (* mixer + lowpass i *)
    r2:=r2 + (r3-r2)*bw;
    r3:=r3 + (r1-r2)*bw2;
    i1:=i1 + (x.Re*ddsi + x.Im*ddsr -i1-i3)*bw;                      (* mixer + lowpass q *)
    i2:=i2 + (i3-i2)*bw;
    i3:=i3 + (i1-i2)*bw2;
  UNTIL i>=b;
  WITH rx^.tapre DO ucr1:=r1; ucr2:=r2; ilr:=r3; ucr2:=VAL(REAL, r2) END;
  WITH rx^.tapim DO ucr1:=i1; ucr2:=i2; ilr:=i3; ucr2:=VAL(REAL, i2) END;
  rx^.phase:=ph;
END iirvarf4;


PROCEDURE mkfir(fg, gain:REAL; flen:CARDINAL; VAR fir:ARRAY OF REAL);
VAR i, m:CARDINAL;
    u, w, l:REAL;
BEGIN
  w:=(2.0*PI)*fg;
  m:=flen DIV 2;
  fir[m]:=1.0;
  FOR i:=1 TO flen DIV 2 DO
    u:=FLOAT(i)*w;
    u:=sin(u)/u;
    IF m+i<flen THEN fir[m+i]:=u END;
    fir[m-i]:=u;
  END;
  u:=2.0/FLOAT(flen-1);
  l:=0.0;
  FOR i:=0 TO flen-1 DO                                      (* hamming *)
    fir[i]:=fir[i]*(0.54+0.46*cos(PI*(FLOAT(i)*u-1.0)));
    l:=l+fir[i];
  END;
  l:=gain/l;
  FOR i:=0 TO flen-1 DO fir[i]:=fir[i]*l END;                (* gain correction *)

--FOR i:=0 TO flen-1 DO WrFixed(fir[i], 4,1); WrStr(" ") END;
--WrStrLn("");WrFixed(fg, 4,1);WrStrLn("");

END mkfir;

PROCEDURE genfir(fg, gain:REAL; tablen:CARDINAL; VAR f:pFIR);
CONST MINTAB=4;
      MINF=0.001;
VAR i:CARDINAL;
    h:REAL;
    on:BOOLEAN;
BEGIN
  on:=(tablen>=MINTAB) & (fg>MINF);
  IF f=NIL THEN
    IF on THEN ALLOCATE(f, SIZE(f^)) END;                        (* new filter *)
  ELSIF NOT on THEN DEALLOCATE(f, SIZE(f^)); f:=NIL END;         (* remove filter *)

  IF (f<>NIL) & ((f^.lastlen<>tablen) OR (f^.fg<>fg)) THEN       (* parameters changed *)
    f^.fg:=fg;
    f^.lastlen:=tablen;
    i:=HIGH(f^.tab)+1;
    IF tablen>i THEN tablen:=i END;
    tablen:=(tablen+3) DIV 4*4-1;
    f^.len:=tablen;
    f^.halfband:=0;
    IF ABS(fg-0.5)<0.001 THEN 
      fg:=0.5;
      f^.halfband:=tablen DIV 2;
    END;
    f^.wp:=0;
    mkfir(fg*0.5, gain, tablen, f^.tab);
(*
FOR i:=0 TO tablen-1 DO WrFixed(f^.tab[i], 4,1); WrStr(" ") END;
WrStrLn("");WrFixed(f^.tab[f^.halfband], 4,1); WrStrLn("");WrInt(f^.halfband,1);WrStrLn("=half");
*)
  END;
END genfir;


PROCEDURE dofir(VAR u:Complex; fir:pFIR; flt:BOOLEAN);
VAR i,j:CARDINAL;
    t:REAL;
BEGIN
  WITH fir^ DO
    re[wp]:=u.Re;
    im[wp]:=u.Im;
    INC(wp);
    IF flt THEN
      j:=0; 
      IF halfband>0 THEN
        i:=(wp+halfband) MOD len;
        t:=tab[halfband];
        u.Re:=re[i]*t;
        u.Im:=im[i]*t;
        i:=wp;
        WHILE i<len DO
          t:=tab[j];
          u.Re:=u.Re+re[i]*t;
          u.Im:=u.Im+im[i]*t;
          INC(j,2);
          INC(i,2);
        END;
        i:=0;
        WHILE i<wp DO
          t:=tab[j];
          u.Re:=u.Re+re[i]*t;
          u.Im:=u.Im+im[i]*t;
          INC(j,2);
          INC(i,2);
        END;
      ELSE
        u.Re:=0.0;
        u.Im:=0.0;
        i:=wp;
        WHILE i<len DO
          t:=tab[j];
          u.Re:=u.Re+re[i]*t;
          u.Im:=u.Im+im[i]*t;
          INC(j);
          INC(i);
        END;
        i:=0;
        WHILE i<wp DO
          t:=tab[j];
          u.Re:=u.Re+re[i]*t;
          u.Im:=u.Im+im[i]*t;
          INC(j);
          INC(i);
        END;
      END;
    END;
    wp:=wp MOD len;
  END;
END dofir;


PROCEDURE ssbiir(VAR tap:SSBTAP; fg, fgq:REAL; VAR u:REAL);
VAR i,j:CARDINAL;
BEGIN
  WITH tap DO
    uc1:=uc1 + (u-uc1)*fg - il;
    uc2:=uc2 + il-il2;
    uc3:=uc3 + il2-uc3*fg;
    il :=il  + (uc1-uc2)*fgq;
    il2:=il2 + (uc2-uc3)*fgq;
    u:=uc3;
  END;
END ssbiir;


PROCEDURE rotvector(VAR v:Complex; rr, ri:REAL);
VAR t:REAL;
BEGIN
  t   :=v.Re*rr - v.Im*ri;
  v.Im:=v.Re*ri + v.Im*rr;
  v.Re:=t;
END rotvector; 


PROCEDURE getsamp(rx:pRX; notfirst:BOOLEAN):INT16;
VAR u, abs:Complex;
    w, af:REAL;
    lev,l:REAL;
BEGIN
  u.Re:=rx^.tapre.ucr2;
  u.Im:=rx^.tapim.ucr2;

(* fine shift rest of full 1khz *)
  IF rx^.fine>0 THEN
    rx^.fracphase:=CAST(CARDINAL, CAST(SET32, rx^.fracphase + rx^.fine) * CAST(SET32, SSBDDSSIZE-1));
    rotvector(u, SSBDDS[rx^.fracphase],
    SSBDDS[CAST(CARDINAL, CAST(SET32, rx^.fracphase + SSBDDSSIZE4) * CAST(SET32,SSBDDSSIZE-1))]);
  END;
(* fine shift rest of full 1khz *)


  IF rx^.modulation=mSSB THEN
(* ssb *)

(* additional IF fir *)
    ssbiir(rx^.ssbre, rx^.ssbfg, rx^.ssbfgq, u.Re);
    ssbiir(rx^.ssbim, rx^.ssbfg, rx^.ssbfgq, u.Im);
(* additional IF fir *)

(* rssi *)
    lev:=1.0 + u.Re*u.Re + u.Im*u.Im;
    l:=lev-rx^.rssi;
    IF l>=0.0 THEN l:=l*0.1 ELSE l:=l*rx^.agcspeed END;                 (* fast up slow down *)
    rx^.rssi:=rx^.rssi + l;
(* rssi *)

(* ssb *)
    
    rx^.bfophase:=CAST(CARDINAL, CAST(SET32, rx^.bfophase + rx^.bfo) * CAST(SET32, SSBDDSSIZE-1));
    af:=(u.Re*SSBDDS[rx^.bfophase]
       - u.Im*SSBDDS[CAST(CARDINAL, CAST(SET32, rx^.bfophase + SSBDDSSIZE4) * CAST(SET32,SSBDDSSIZE-1))])
       *25000.0/sqrt(rx^.rssi);
(* ssb *)

  ELSE                                           (* AM FM *)

    IF rx^.modulation=mSCAN THEN

(* scan squelch *)
      lev:=u.Re*u.Re + u.Im*u.Im;
      IF notfirst THEN rx^.sqsum:=rx^.sqsum + ABS(rx^.lastlev-lev);    (* sum noise *)
      ELSE
        rx^.sqsum:=0.0;
        rx^.rssi:=0.0;
      END;
      rx^.rssi:=rx^.rssi + lev;                                        (* sum levels *)
      rx^.lastlev:=lev;
(* scan squelch *)

      af:=sqrt(lev)*0.01;
    ELSE

      IF rx^.fir<>NIL THEN dofir(u, rx^.fir, TRUE) END;

(* rssi *)
      lev:=1.0 + u.Re*u.Re + u.Im*u.Im;
      l:=lev-rx^.rssi;
      rx^.rssi:=rx^.rssi + l*0.001;
(* rssi *)

(* complex to phase *)
      abs.Re:=ABS(u.Re);
      abs.Im:=ABS(u.Im);
      IF abs.Im>abs.Re THEN
        IF abs.Im>0.0 THEN w:=abs.Re/abs.Im ELSE w:=0.0 END;
        w:=pi/2.0 - (w*1.055 - w*w*0.267);                    (* arctan *)
      ELSE
        IF abs.Re>0.0 THEN w:=abs.Im/abs.Re ELSE w:=0.0 END;
        w:=w*1.055 - w*w*0.267;
      END;
      IF u.Re<0.0 THEN w:=pi-w END;
      IF u.Im<0.0 THEN w:=-w END;
(* complex to phase *)

(* phase highpass make FM *)
      af:=w-rx^.w1;
      rx^.w1:=w;
      IF af> pi THEN af:=af - pi*2.0 END;
      IF af<-pi THEN af:=af + pi*2.0 END;
(* phase highpass make FM *)

      IF rx^.modulation=mAM THEN

(* am squelch *)
        IF rx^.squelch THEN
          IF notfirst THEN rx^.sqsum:=rx^.sqsum + ABS(rx^.a1-af) END;
          rx^.a1:=af;
        END;
(* am squelch *)

(* am demod *)
        lev:=sqrt(lev);                         (* amplitude *)
        af:=lev-rx^.lastlev;                    (* - median aplitude *) 
        rx^.lastlev:=rx^.lastlev + af*0.001;
        af:=af/rx^.lastlev*32000.0;             (* agc *)
(* am demod *)

      ELSE

(* fm squelch *)
        IF rx^.squelch THEN
          IF notfirst THEN rx^.sqsum:=rx^.sqsum + ABS(rx^.lastlev-lev)/(rx^.lastlev+lev) END;
          rx^.lastlev:=lev;
        END;
(* fm squelch *)

        af:=af*(25000.0/pi);

      END;
    END;
  END;
  IF af>32000.0 THEN af:=32000.0 ELSIF af<-32000.0 THEN af:=-32000.0 END;

  RETURN VAL(INT16, af);
END getsamp;

PROCEDURE getiq(rx:pRX; VAR i,q:INT16; dofilt:BOOLEAN);
CONST M=32767.0;
VAR u:Complex;
BEGIN
  u.Re:=rx^.tapre.ucr2*iqlevel;
  u.Im:=rx^.tapim.ucr2*iqlevel;
  IF rx^.fir<>NIL THEN dofir(u, rx^.fir, dofilt) END;
  IF dofilt THEN
    IF u.Re>M THEN u.Re:=M ELSIF u.Re<-M THEN u.Re:=-M END;
    IF u.Im>M THEN u.Im:=M ELSIF u.Im<-M THEN u.Im:=-M END;
    i:=VAL(INT16, u.Re);
    q:=VAL(INT16, u.Im);
  END;
END getiq;


PROCEDURE getsdr(samps:CARDINAL; rx-:ARRAY OF pRX):INTEGER;
CONST FINESTEP=1024;

VAR s,r,a,b,ws,wssb, inwords, insamps, inbytes:CARDINAL;
    u,res:INTEGER;
    allsamp:REAL;
BEGIN
  IF reconnect & (fd<0) THEN
    IF isfile THEN fd:=OpenRead(url)
    ELSE
      usleep(1000000);
      fd:=connecttob(url,port);
    END;
  END;
  IF fd>=0 THEN 
    IF samps>IQBUF DIV 2-1 THEN samps:=IQBUF DIV 2-1 END;
    allsamp:=sampfrac + onesamp*FLOAT(samps);
    insamps:=TRUNC(allsamp);
    inwords:=insamps*2;
    inbytes:=insamps*bytespersamp;

    IF isfile THEN
      a:=0; 
      REPEAT
        res:=read(fd, ADR(iqbuf.b[a]), inbytes-a);
        IF res<=0 THEN
          Close(fd);
          fd:=-1;
          RETURN -1

        END;
        INC(a, res);
      UNTIL a>=inbytes;

    ELSIF readsockb(fd, iqbuf, inbytes)<0 THEN  (* connect lost *)
      Close(fd);
      fd:=-1;
      RETURN -1

    END;
    IF debfd>=0 THEN WrBin(debfd, iqbuf, inbytes) END;

    IF bytespersamp=2 THEN
      FOR a:=insamps-1 TO 0 BY -1 DO
        iqbuf.i32[a].r:=VAL(INTEGER, iqbuf.u8[a].Re)-127;
        iqbuf.i32[a].i:=VAL(INTEGER, iqbuf.u8[a].Im)-127;
      END;     
    ELSIF bytespersamp=4 THEN
      FOR a:=insamps-1 TO 0 BY -1 DO
        iqbuf.i32[a].r:=VAL(INTEGER, iqbuf.i16[a].Re);
        iqbuf.i32[a].i:=VAL(INTEGER, iqbuf.i16[a].Im);
      END;
    END;

    a:=0;
    FOR s:=0 TO samps-1 DO
      r:=0;
      b:=TRUNC(sampfrac + onesamp*FLOAT(s+1));
      WHILE rx[r]<>NIL DO
        WITH rx[r]^ DO
          ws:=width;
          IF modulation=mSSB THEN
            wssb:=width*5 DIV 4;
            ws:=3000;
            WHILE ws<wssb DO ws:=ws*2 END;
--WrInt(ws, 0);WrStrLn(" ws");
            ssbfg:=FLOAT(width)/FLOAT(audiohz);
            ssbfgq:=ssbfg*ssbfg*2.0;
            bfo:=VAL(CARDINAL, ABS(maxafc))*SSBDDSSIZE DIV audiohz;
            IF maxafc>0 THEN bfo:=SSBDDSSIZE-bfo END;
--WrInt(bfo, 0);WrStrLn(" bfo");
            fine:=dffrac*SSBDDSSIZE DIV audiohz;

            IF agc>0 THEN agcspeed:=0.2/FLOAT(agc) ELSE agcspeed:=0.00025 END;
          ELSE fine:=0; agcspeed:=0.001 END;

          IF bytespersamp=8 THEN iirvarf4(rx[r], a, b, FLOAT(width)*varbwk);
          ELSE
            IF rtlhz<2048000 THEN ws:=ws*2 END; 
            IF ws=3000 THEN iir512(rx[r], a, b);
            ELSIF ws=6000 THEN iir256(rx[r], a, b);
            ELSIF ws=12000 THEN iir128(rx[r], a, b);
            ELSIF ws=24000 THEN iir64(rx[r], a, b);
            ELSIF ws=48000 THEN iir32(rx[r], a, b);
            ELSIF ws=96000 THEN iir16(rx[r], a, b);
            ELSIF ws=192000 THEN iir8(rx[r], a, b);
            ELSE iirvar(rx[r], a, b, FLOAT(ws)*(1.0/(11500.0*128.0))) END;
          END;
          IF modulation=mIQ THEN
            IF qsamples<>NIL THEN getiq(rx[r], samples^[s], qsamples^[s], TRUE) (* full rate iq *)
            ELSE getiq(rx[r], samples^[s], samples^[s+1], NOT ODD(s)) END; (* half rate iq *)
          ELSE
            u:=getsamp(rx[r], s>0);
            samples^[s]:=u;
            IF afcrun THEN median:=median + u END;                             (* afc *)  
          END;
        END;
        INC(r);
      END;
      a:=b;
    END;
    sampfrac:=allsamp-FLOAT(insamps);

--AFC
    r:=0;
    WHILE rx[r]<>NIL DO
      WITH rx[r]^ DO
        afcrun:=FALSE; 
        IF (modulation=mFM) & (maxafc>0) THEN
          IF afcrssi<rssi THEN afcrssi:=rssi ELSE afcrssi:=afcrssi*AFCRSSISPEED END; (* peak level rssi *)
          IF afcrssi*AFCRSSICLAMP<rssi THEN                                     (* near peak rssi so let afc work *)
            median:=median - afckhz*VAL(INTEGER, samps)*AFCRECENTER DIV maxafc; (* weak pull to middle *)
            IF median>afcspeed THEN
              INC(afckhz);
              IF afckhz>maxafc THEN afckhz:=maxafc END;
              median:=0;
            ELSIF median<-afcspeed THEN
              DEC(afckhz);
              IF afckhz<-maxafc THEN afckhz:=-maxafc END;
              median:=0;
            END;
            afcrun:=TRUE;
          END;                                                   (* low rssi so freeze afc *)
        END;
--WrInt(ORD(afcrun),1);
      END;
      INC(r);
    END;
--AFC

    RETURN samps
    
  ELSE RETURN -1 END;
END getsdr;

PROCEDURE setparm(num, value:CARDINAL);
VAR res:INTEGER;
    tbuf:ARRAY[0..4] OF CHAR;
BEGIN
  tbuf[0]:=CHR(num);
  tbuf[1]:=CHR(value DIV 1000000H);
  tbuf[2]:=CHR(value DIV 10000H MOD 100H);
  tbuf[3]:=CHR(value DIV 100H MOD 100H);
  tbuf[4]:=CHR(value MOD 100H);
  res:=sendsock(fd, tbuf, 5);
END setparm;


PROCEDURE startsdr(ip-,tport-:ARRAY OF CHAR; inhz, outhz:CARDINAL; reconn:BOOLEAN; format:CARDINAL):BOOLEAN;
BEGIN
  Assign(url, ip);
  Assign(port, tport);
  reconnect:=reconn;
  bytespersamp:=2;
  iqlevel:=IQLEVEL;
  IF format=2 THEN
    bytespersamp:=4;
    iqlevel:=IQLEVEL/256.0;
  END;
  IF inhz>0 THEN rtlhz:=inhz END;
  IF format=4 THEN
    IF (rtlhz<1000000) OR (rtlhz>3000000) THEN RETURN FALSE END;
    bytespersamp:=8;
    iqlevel:=IQLEVEL*128.0;
    varbwk:=2048000.0/(11500.0*128.0)/FLOAT(rtlhz);
  ELSIF (rtlhz<>1024000) & ((rtlhz<2048000) OR (rtlhz>2500000)) THEN RETURN FALSE END;

  IF outhz>0 THEN audiohz:=outhz END;
  onesamp:=FLOAT(rtlhz)/FLOAT(audiohz);                            (* sample reduction *)
  sampfrac:=0.0;
  isfile:=(port[0]="0") & (port[1]=0C);
  IF fd<0 THEN
    IF isfile THEN fd:=OpenRead(url) ELSE fd:=connecttob(url,port) END;
  END;
  IF fd>=0 THEN
    IF NOT isfile & (format=1) THEN setparm(2, rtlhz) END;
    IF inhz>=2048000 THEN initdds(2048) ELSE initdds(1024) END;
    initssbdds(SSBDDS);
  END;
  RETURN fd>=0
END startsdr;


BEGIN
  fd:=-1;
  debfd:=-1;
  reconnect:=FALSE;
  rtlhz:=2048000;
  audiohz:=16000;
  afcspeed:=AFCSPEED;
END sdr.