{$A+ word align data}
{$G+ 286 instructions}
{$R- range checking}
{$S- stack checking}
{$Q- overflow checking}

uses crt,graph;
{$L EGAVGA} { Чтобы не таскать за собой egavga.bgi     }
            { Для получения файла egavga.obj выполнить }
            { binobj egavga.bgi egavga egavgabgi       }


{$define correct} {коppектиpовать отpицательные после деления на 2  }
                  {иначе всегда будет округляться в меньшую сторону,}
                  {а лучше - чтобы к нулю - повышается устойчивость }

procedure egavgabgi; external;

   const n=2*256;  {число элементов в массиве = 2*число отсчетов}
         scale=2;  {Масштаб для регулировки уровня входного сигнала.
                    Стоит подбирать или сотворить АРУ, т.к. мало бит
                    при расчетах}
TYPE
   gldarray =^tgldarray;
   tgldarray = ARRAY [0..N] OF real; {по очереди вещественные и мнимые
                                      составляющие}
   igldarray =^tigldarray;
   tigldarray = ARRAY [0..N] OF integer; { то же, но в целых }
 TYPE double = real;
var    cost,sint:tigldarray;


PROCEDURE four1( data: gldarray; nn,isign: integer);
  {БПФ в real; списан с какой-то книжки}
(* Programs using routine FOUR1 must define type*)
(*in the calling routine, where nn2=nn+nn. *)
VAR
   ii,jj,n,mmax,m,j,istep,i: integer;
   wtemp,wr,wpr,wpi,wi,theta: real;
   tempr,tempi: real;
BEGIN
   n := 2*nn;
   j := 1;
   FOR ii := 1 TO nn DO BEGIN
      i := 2*ii-1;
      IF (j > i) THEN BEGIN
         tempr := Data^[j];
         tempi := Data^[j+1];
         Data^[j] := Data^[i];
         Data^[j+1] := Data^[i+1];
         Data^[i] := tempr;
         Data^[i+1] := tempi
      END;
      m := n DIV 2;
      WHILE ((m >= 2) AND (j > m))  DO BEGIN
         j := j-m;
         m := m DIV 2
      END;
      j := j+m
   END;
   mmax := 2;
   WHILE (n > mmax) DO BEGIN
      istep := 2*mmax;
      theta := 6.28318530717959/(isign*mmax);
      wpr := -2.0*sqr(sin(0.5*theta));
      wpi := sin(theta);
      wr := 1.0;
      wi := 0.0;
      FOR ii := 1 TO (mmax DIV 2) DO BEGIN
         m := 2*ii-1;
         FOR jj := 0 TO ((n-m) DIV istep) DO BEGIN
            i := m + jj*istep;
            j := i+mmax;
            tempr := trunc(wr*Data^[j]-wi*Data^[j+1]);
            tempi := trunc(wr*Data^[j+1]+wi*Data^[j]);
            Data^[j] := Data^[i]-tempr;
            Data^[j+1] := Data^[i+1]-tempi;
            Data^[i] := Data^[i]+tempr;
            Data^[i+1] := Data^[i+1]+tempi
         END;
         wtemp := wr;
         wr := wr*wpr-wi*wpi+wr;
         wi := wi*wpr+wtemp*wpi+wi
      END;
      mmax := istep
   END
END;

function LongMul(X, Y: Integer): Longint; inline($5A/$58/$f7/$EA);
function LongDiv(X: Longint; Y: Integer): Integer; inline($59/$58/$5A/$F7/$F9);

PROCEDURE fourtrunc( data: igldarray; nn,isign: integer);
(* Programs using routine FOUR1 must define type*)
(*in the calling routine, where nn2=nn+nn. *)
VAR
   ii,jj,n,mmax,m,j,istep,i: integer;
   wr,wi,fi:integer;
   tempr,tempi: integer;
label m0,m1,m2,m3,m4;
BEGIN
   n := 2*nn;
   j := 1;
   FOR ii := 1 TO nn DO BEGIN
      i := 2*ii-1;
      IF (j > i) THEN BEGIN
         tempr := Data^[j];
         tempi := Data^[j+1];
         Data^[j] := Data^[i];
         Data^[j+1] := Data^[i+1];
         Data^[i] := tempr;
         Data^[i+1] := tempi
      END;
      m := n shr 1;
      WHILE ((m >= 2) AND (j > m))  DO BEGIN
         j := j-m;
         m := m DIV 2
      END;
      j := j+m
   END;
   mmax := 2;
   WHILE (n > mmax) DO BEGIN
      istep := mmax shl 1;
{}    fi := 0;
      FOR ii := 1 TO mmax shr 1 DO BEGIN
         wr := cost[fi];
         wi := sint[fi];
         m := ii shl 1-1;
         i:=m;
          asm
{           FOR jj := 0 TO ((n-m) DIV istep) DO BEGIN}
            mov  ax,n
            sub  ax,m
            xor  dx,dx
            div  istep
            inc  ax
            mov  jj,ax { беpем ..-jj }
m0:
              push ds
              mov  ax,word ptr data+2 {seg data^}
              mov  ds,ax

              mov  si,i
              mov  di,si
{             j := i+mmax;}
              add  di,mmax
              add  si,si
              add  si,word ptr data {offset data^}
              add  di,di
              add  di,word ptr data {offset data^}

{            tempr :=
longdiv((longMul(wr,Data^[j])-longmul(wi,Data^[j+1])),128);{}
              mov  ax,[ds:di]
              imul wr
              sal  ax,1
              rcl  dx,1
              mov  bl,ah
              mov  bh,dl
              mov  ax,[ds:di+2]
              imul wi
              sal  ax,1
              rcl  dx,1
              mov  al,ah
              mov  ah,dl
              sub  bx,ax
              mov  tempr,bx

{            tempi :=
longdiv((longmul(wr,Data^[j+1])+longmul(wi,Data^[j])),128);{}
              mov  ax,[ds:di]
              imul wi
              sal  ax,1
              rcl  dx,1
              mov  cl,ah
              mov  ch,dl
              mov  ax,[ds:di+2]
              imul wr
              sal  ax,1
              rcl  dx,1
              mov  al,ah
              mov  ah,dl
              add  cx,ax
              mov  tempi,cx

{            Data^[j] := (Data^[i]-tempr) div 2;}
              mov  ax,[ds:si]
              sub  ax,bx {tempr}
              sar  ax,1
              {$ifdef correct}
              jns  m1 { пpи отpицательных значениях окpугление идет не так, как
с div }
              adc  ax,0
              {$endif}
m1:           mov  [ds:di],ax

{            Data^[j+1] := (Data^[i+1]-tempi) div 2;}
              mov  ax,[ds:si+2]
              sub  ax,cx {tempi}
              sar  ax,1
              {$ifdef correct}
              jns  m2 { пpи отpицательных значениях окpугление идет не так, как
с div }
              adc  ax,0
              {$endif}
m2:           mov  [ds:di+2],ax

{            Data^[i] := (Data^[i]+tempr) div 2;}
              mov  ax,[ds:si]
              add  ax,bx {tempr}
              sar  ax,1
              {$ifdef correct}
              jns  m3 { пpи отpицательных значениях окpугление идет не так, как
с div }
              adc  ax,0
              {$endif}
m3:           mov  [ds:si],ax

{            Data^[i+1] := (Data^[i+1]+tempi) div 2;}
              mov  ax,[ds:si+2]
              add  ax,cx {tempi}
              sar  ax,1
              {$ifdef correct}
              jns  m4 { пpи отpицательных значениях окpугление идет не так, как
с div }
              adc  ax,0
              {$endif}
m4:           mov  [ds:si+2],ax

              pop ds
{             inc(i,istep);}
              mov  ax,istep
              add  i,ax
{           END;}
            dec  jj
            jnz  m0
          end;
{}       fi := fi + nn*2 div mmax;
      END;
      mmax := istep
   END
END;

var t:tgldarray;
    q,ttrunc,t1:tigldarray;
    i,j:word;
    time:longint; {таймер для подсчета времени работы}
    number:word;
  var BIOSTicks : longint absolute $40:$6C; {таймер доса}
      f:file; {input file; stereo 8bit}
      b:array[1..n] of byte;
      tmp:integer;
var
 grDriver: Integer;
 grMode: Integer;
 ErrCode: Integer;

Procedure PlotPixel( x, y, i: word ); {вывод точки спектра на экран соотв.
                                       цветом}
const ColNum = 9; { число цветов }
      OffTo = 3; { до такого i не отобpажать }
const colors:array[0..ColNum] of byte =
     ( Cyan,Black,Red,Red,Red,Brown,LightRed,LightRed,LightRed,Yellow );
begin
  if i>OffTo then i := i-OffTo else i := 0;
  if i>ColNum then i := ColNum;
  PutPixel( x, GetMaxY - y, Colors[i] );
end;

function energy(x,y,freq:integer):word; { энеpгия от косинусной и синусной
                                          амплитуд и частоты }
var i,w:word;
begin
  {Если не нужна большая точность, то вместо
   вычисления квадратного корня
   w:=round(sqrt(sqr(x)+sqr(y))
   быстрее сделать так}
  x:=abs(x);
  y:=abs(y);
  if      x<y shr 1 then w:=x
  else if y<x shr 1 then w:=y
                    else w:=(x+y) shr 1;
  if w<128 then energy:=w else energy:=127;
  {это бы все тоже лучше на асм}
end;

var r,l:array[1..n div 4] of byte;

begin
  grDriver := Detect;
  RegisterBgiDriver(@egavgabgi);
  InitGraph(grDriver, grMode,' ');
  ErrCode := GraphResult;
  if ErrCode <> grOk then begin
    Writeln('Graphics error:', GraphErrorMsg(ErrCode));
    halt(254);
  end;

  { для подсчета времени работы процедур }
  time:=0; {обнулим таймер}
  number:=0; {кол-во квантов}

  assign(f,'handsoff.wav');reset(f,1); {файл, из которого читаем;
                                        стерео 8 бит}
  for i:=0 to N do begin { заполним таблицы синусов и косинусов }
    cost[i]:=round(128*cos(2*pi*i/n));
    sint[i]:=round(128*sin(2*pi*i/n));
  end;

  j:=0; {позиция столбца изображения спектра на экране}

 repeat
{$i-}  blockread(f,b,n); {$i+}
  if ioresult<>0 then begin
    reset(f,1);
    setcolor(black);
    line(j-1,0,j-1,getmaxy);
  end;

{ для подсчета времени работы процедур }
  time :=time-biosticks; { включить таймер}

  for i:=1 to N do begin
    ttrunc[n+1-i]:=b[i]*Scale; {Почему-то в обратном порядке.
                                Абсолютно не помню, зачем}
  end;

  fourtrunc(@ttrunc,n div 2,1);
{  four1(@t,n div 2,1);}

  for i:=1 to N do b[i]:=0;

  for i:=1 to n div 4 - 1 do begin
    {Коэфф. при exp(i*k*t) приводятся к коэф. при sin и cos,
     разделяем левый и правый каналы}
    tmp := ttrunc[i*2+2]+ttrunc[n+2-i*2];          {косинусы 1 канала}
    ttrunc[i*2+2]:=ttrunc[i*2+2]-ttrunc[n+2-i*2];  {синусы 1 канала}
    ttrunc[n+2-i*2]:=tmp;                          {косинусы 1 канала}

    tmp := ttrunc[i*2+1]-ttrunc[n+1-i*2];          {синусы 2 канала}
    ttrunc[i*2+1]:=ttrunc[i*2+1]+ttrunc[n+1-i*2];  {косинусы 2 канала}
    ttrunc[n+1-i*2]:=tmp;                          {синусы 2 канала}

    {я уже подзабыл, в каком порядке раскладываются sin и cos,
     но вроде по порядку}

    r[i]:=energy(ttrunc[i*2+2],ttrunc[i*2+1],i);
    l[i]:=energy(ttrunc[n+2-i*2],ttrunc[n+1-i*2],i);
    {энергетические функции для 2-х каналов}

  end;

  time :=time+biosticks; {выключить таймер}
  for i:=1 to n div 4 do begin
    plotpixel(j,i+280,l[i]);

    plotpixel(j,i+60,r[i]);
  end;

  setcolor(white);

  inc(j); if j>getmaxx then j:=0;

  inc(number);

  until keypressed;
  readkey; closegraph;

(* {Если есть желание смотреть результат в
    сравнении с точным - раскомментируй}
{$i-}  blockread(f,b,n); {$i+}
  for i:=1 to N do begin
    ttrunc[n+1-i]:=b[i]*Scale;
    t[n+1-i]:=ttrunc[n+1-i];
  end;

  four1(@t,n div 2,1);
  fourtrunc(@ttrunc,n div 2,1);

  writeln;
  writeln;
  for i:=1 to N do begin
    t1[i]:=round(t[i]/256);
    write(t1[i]:4,'=',ttrunc[i]:4,' ');
  end; {чтобы сравнить с точным результатом} (* *)
  writeln; writeln(number,' times,
',time/18.2:6:2,'sec,',number/time*18.2:3:0,' times/sec.');{}
end.

