DEMO.DESIGN Frequently Asked Questions |  | Интернет: ЛЭИВО

оглавление | demo party в ex-СССР | infused bytes e-mag | новости от ib/news | другие проекты | письмо | | win koi lat

    следующий фpагмент (2) From : Andrew Zabolotny                    2:5030/84.5     28 Jan 96  21:18:50
     To   : Anatoly Davidov                     5075/18.2       29 Jan 96  20:32:38
     Subj : Real-time spectrum analiser..
    --------------------------------------------------------------------------------
    Hello Anatoly!

    Saturday January 27 1996 19:40, Anatoly Davidov wrote to Andrey Zabolotny:

     AD>  Вот имеется у меня твой Player, а в нем есть приятная такая весчь - subj!
     AD> Который интересует меня чисто как subj. Ведь вроде бы использовался
     AD> алгоритм быстрого преобразования Фурье, не так ли ? Вот недавно сделал
     AD> АЦП, и очень хотелось бы поэксперементировать над саундом. Вот. А не мог
     AD> бы ты случайно поделится этим алгоритмом, или может исходнички какие не
     AD> жалко..? Или может ты использовал какие статьи из журналов, (говорят в
     AD> каком-то пробегало подобное), то оочень тебя прошу, тсказать, шепнуть
     AD> номерок. А? Пожалуйста!

    Hу вообще-то player у меня давно уже похеpен. Hо чтой-то такое недавно делал.
    Ща поищу... да, вот он - ниже.

    Главная идея такая: для каждой частоты котоpую тебе нужно отобpазить сочиняешь
    таблицу sin(Freq) и cos(Freq). То есть как если бы мы оцифpовали сигнал с
    опpеделенной амплитудой и частотой Freq на твоей частоте дискpетизации.
    Длина таблицы должна быть больше чем одна полуволна сигнала Freq. То есть если
    ты делаешь несколько таблиц для нескольких частот то длина каждой (все таблицы
    pавны в длину!) беpешь длину полуволны для низшей частоты.
    Далее для опpеделения спектpа кусочка сигнала (с длиной pавной длине таблиц)
    для каждой частоты анализатоpа спектpа считаешь сумму пpоизведений:
    A := Сумма(Signal[I] * Sin[I]) для I=(1..BlockLen)
    B := Сумма(Signal[I] * Cos[I]) для I=(1..BlockLen)
    После этого амплитуда данной частоты есть
    Sqrt(A^2 + B^2).
    Чтобы не вычислять коpень можешь делать пpосто A^2 + B^2.
    Hу и как-то их масштабиpуешь. Потом pисуешь палку соотв. длины и усе.

    ------------------------------------------------
    uses crt,strOp,fastFloat,miscUtil,kbConst,vPack,sbIO;

    const bPerM  = 63;
          fBands = 79;
          fLo    = 200;
          fHi    = 4000;
          fDiscr = 8000;

    var   SinCos : array[0..(fBands + 1) * (bPerM + 1) * 2 - 1] of ShortInt;

    Procedure FourierInit;
    var I,J : Integer;
        F,T : Real;
    begin
     For I := 0 to fBands do
         begin
          F := exp(ln(fLo) + (ln(fHi) - ln(fLo)) * I / succ(fBands));
          For J := 0 to bPerM do
              begin
               T := 2 * Pi * J * (F / fDiscr);
               SinCos[(J * succ(fBands) + I) * 2 + 0] := round(System.Sin(T) * 64);
               SinCos[(J * succ(fBands) + I) * 2 + 1] := round(System.Cos(T) * 64);
              end;
         end;
    end;

    Procedure ShowAmplitude(var Buffer : tArrOfWord);
    var I,J,H,A : Integer;
    begin
     For I := 0 to fBands do
         begin
          H := iSqrt(lSqr(Buffer[I * 2]) + lSqr(Buffer[I * 2 + 1])) div 2000;
          For J :=0 to 21 do
              begin
               A := 160 * 24 + I * 2 - J * 160;
               if J < H then mem[$B800:A] := 22
                        else mem[$B800:A] := 32;
               case J of
                00..11 : mem[$B800:A+1] := $0A;
                12..18 : mem[$B800:A+1] := $0E;
                19..21 : mem[$B800:A+1] := $0C;
               end;
              end;
         end;
    end;

    Procedure FourierProceed(var Data; Size : Word);
    var SinCosAmp : array[0..fBands * 2 + 1] of Integer;
    begin
     asm            cld
                    mov     ax,ss
                    mov     es,ax
                    lea     di,SinCosAmp
                    mov     cx,(fBands + 1) * 2
                    xor     ax,ax
                    rep     stosw
                    push    bp
                    les     si,Data
                    lea     bp,SinCosAmp
                    mov     bx,offset SinCos

                    mov     dl,bPerM + 1

    @@nextSrcB:     SegES   lodsb
                    xor     al,80h
                    mov     ch,al
                    mov     cl,(fBands + 1) * 2

    @@doFourier:    mov     al,ch
                    imul    byte ptr [bx]
                    add     word ptr [bp],ax
                    jno @@@
                    nop
                    @@@:
                    inc     bp
                    inc     bp
                    inc     bx
                    dec     cl
                    jnz     @@doFourier
                    sub     bp,(fBands + 1) * 2 * 2

                    dec     dl
                    jne     @@nextSrcB

                    pop     bp
     end;
     ShowAmplitude(pArrOfWord(@SinCosAmp)^);
    end;

    {SoundBlaster stuff}

    var outF : File;
        S    : pSample;
        mid  : boolean;

    Procedure ShowMarker;
    var Q : String;
        I : Word;
    begin
     Q := Strg('-',64);
     I := luDiv(luMul(sbGetPos, 64), S^.Size) + 1;
     if I > 64 then I := 64;
     Q[I] := '$';
     GoToXY(1,3);
     TextAttr := $1E; Write('processing (');
     if Mid then TextAttr := $1C else TextAttr := $1F;
     Mid := False; Write(Q);
     TextAttr := $1E; Write(')'); ClrEOL; Writeln;
    end;

    Procedure printMid(Buffer : Pointer; Size : Word); far;
    begin
     asm push es; push cx; push si; push di; end;
     mid := TRUE;
     asm pop di; pop si; pop cx; pop es; end;
    end;

    Procedure flushBuf(Buffer : Pointer; Size : Word); far;
    begin
     asm push es; push cx; push si; push di; end;
     blockWrite(outF, Buffer^, Size);
     mid := TRUE;
     asm pop di; pop si; pop cx; pop es; end;
    end;

    Procedure testPlay;
    begin
     if not loadSample(S, 'record.raw') then Halt;
     sbSamplingRate := 11000;
     sbSetBuffer(S);
     sbStart(opPlayPCM8 or opLoop, printMid);
     repeat
      ShowMarker;
     until keypressed;
     deallocateBuffer(S);
    end;

    Procedure testRecord;
    begin
     assign(outF, 'record.raw'); Rewrite(outF, 1);
     if not allocateBuffer(S, 32768) then Halt;
     sbSamplingRate := 11000;
     sbSetBuffer(S);
     sbStart(opRecordADPCM2 or opLoop, flushBuf);
     repeat
      ShowMarker;
     until keypressed;
     sbStop;
     deallocateBuffer(S);
     close(outF);
    end;

    Procedure testFourier;
    var ii : Word;
    begin
     FourierInit;
     if not allocateBuffer(S, 32768) then Halt;
     sbSamplingRate := fDiscr;
     sbSetBuffer(S);
     sbStart(opRecordADPCM2 or opLoop, printMid);
     repeat
      ShowMarker;
      ii := sbGetPos - bPerM;
      if (ii < S^.Size) and (ii + bPerM < S^.Size)
         then FourierProceed(S^.Data[ii], 1);
      Delay(50);
     until keypressed;
     sbStop;
     deallocateBuffer(S);
     close(outF);
    end;

    var ii   : Word;

    begin
     TextAttr := $07;
     ClrScr;
     if sbDetect
        then begin
              Writeln('Soundblaster detected at port ', Hex4(sbPort), 'h');
              Writeln('IRQ ',sbIRQ,', DMA ',sbDMA,', DSP
    v',hi(sbVERSION),'.',lo(sbVERSION))
             end
        else begin Writeln('Soundblaster not detected'); Halt; end;
     sbInit;
    { TestRecord;
     While keyPressed do ReadKey;
     TestPlay;}
     TestFourier;
     While keyPressed do ReadKey;
     sbDone;

    (*
     assign(outf,'test.vp'); rewrite(outf,1);
     if not loadSample(S, 'test.raw') then Halt;
     if not vpInit then Halt;
    { vpProcess(S^.Data, S^.Size);}
     FourierInit;

     ii := 0;
     repeat
      FourierProceed(S^.Data[ii], S^.Size);
      inc(ii,succ(bPerM));
      if ii + succ(bPerM) >= S^.Size then II := 0;
      Delay(100);
     until keypressed;
    *)
    end.

    следующий фpагмент (3)|пpедыдущий фpагмент (1)This is from _Numerical Recipes_, adapted from FORTRAN to C.  Just
    for the heck of it, I've also attached the original FORTRAN, since I
    used it as a check against my C coding.  For an explanation of the
    theory, the code, and how to use it, see Press, _et_al_, _Numerical
    Recipes_ (ISBN 0-521-30811-9), pp. 449-453.

    							spl

    		  -------------cut here-------------
    #include <stdio.h>
    #include <math.h>

    /*
     * N dimensional DFT, filched from Numerical Recipes and translated
     * by hand from the original FORTRAN.
     */

    #define PI 3.14159265358979323846
    #define TWO_PI ( PI * 2.0 )

    void ndfft( double *data, int *nn, int ndim, int isign )

    {

        int ntot = 1;
        int nprev = 1;
        int idim;
        double i2pi = isign * TWO_PI;

        for ( idim = 0; idim < ndim; idim++ )
    	ntot *= *( nn + idim );

        for ( idim = 0; idim < ndim; idim++ ) {

    	int n = *( nn + idim );
    	int nrem = ntot / ( n * nprev );
    	int ip1 = 2 * nprev;
    	int ip2 = ip1 * n;
    	int ip3 = ip2 * nrem;
    	int i2rev = 0;
    	int i2;
    	int ifp1;

    	/*
    	 * Bit reversal stuff.
    	 */

    	for ( i2 = 0; i2 < ip2; i2 += ip1 ) {

    	    int ibit;

    	    if ( i2 < i2rev ) {

    		int i1;

    		for ( i1 = i2; i1 < i2 + ip1; i1 += 2 ) {

    		    register int i3;

    		    for ( i3 = i1; i3 < ip3; i3 += ip2 ) {

    			register int i3rev = i2rev + i3 - i2;
    			register double tempr = *( data + i3 );
    			register double tempi = *( data + i3 + 1 );

    			*( data + i3 ) = *( data + i3rev );
    			*( data + i3 + 1 ) = *( data + i3rev + 1 );
    			*( data + i3rev ) = tempr;
    			*( data + i3rev + 1 ) = tempi;

    		    }

    		}

    	    }

    	    ibit = ip2 / 2;
    	    while ( ( ibit > ip1 ) && ( i2rev > ibit - 1 ) ) {

    	        i2rev -= ibit;
    		ibit /= 2;

    	    }

    	    i2rev += ibit;

    	}

    	/*
    	 * Danielson-Lanczos stuff.
    	 */

    	ifp1 = ip1;
    	while ( ifp1 < ip2 ) {

    	    int ifp2 = 2 * ifp1;
    	    double theta = i2pi / ( ( double ) ifp2 / ip1 );
    	    register double wpr;
    	    register double wpi;
    	    register double wr = 1.0;
    	    register double wi = 0.0;
    	    int i3;

    	    wpr = sin( 0.5 * theta );
    	    wpr *= wpr * -2.0;
    	    wpi = sin( theta );

    	    for ( i3 = 0; i3 < ifp1; i3 += ip1 ) {

    		register int i1;
    		register double wtemp;

    		for ( i1 = i3; i1 < i3 + ip1; i1 += 2 ) {

    		    register int i2;

    		    for ( i2 = i1; i2 < ip3; i2 += ifp2 ) {

    			register int i21 = i2 + 1;
    			register int k2 = i2 + ifp1;
    			register int k21 = k2 + 1;
    			register double tempr = ( wr * *( data + k2 ) ) -
    			    ( wi * *( data + k21 ) );
    			register double tempi = ( wr * *( data + k21 ) ) +
    			    ( wi * *( data + k2 ) );

    			*( data + k2 ) = *( data + i2 ) - tempr;
    			*( data + k21 ) = *( data + i21 ) - tempi;

    			*( data + i2 ) += tempr;
    			*( data + i21 ) += tempi;

    		    }

    		}

    		wtemp = wr;
    		wr += ( wr * wpr ) - ( wi * wpi );
    		wi += ( wi * wpr ) + ( wtemp * wpi );

    	    }
    	    ifp1 = ifp2;

    	}
    	nprev *= n;

        }

    }
    			----- and here ------
          subroutine fourn( data, nn, ndim, isign )

          implicit double precision ( a-h, o-z )
          double precision data(*)
          integer nn(ndim)
          integer ndim
          integer isign

          ntot = 1

          do idim = 1, ndim

             ntot = ntot * nn(idim)

          enddo

          nprev = 1

          do idim = 1, ndim

             n = nn(idim)

             nrem = ntot / ( n * nprev )

             ip1 = 2 * nprev
             ip2 = ip1 * n
             ip3 = ip2 * nrem

             i2rev = 1

             do i2 = 1, ip2, ip1

                if ( i2 .lt. i2rev ) then

                   do i1 = i2, i2 + ip1 - 2, 2

                      do i3 = i1, ip3, ip2

                         i3rev = i2rev + i3 - i2

                         tempr = data(i3)
                         tempi = data(i3 + 1)

                         data(i3) = data(i3rev)
                         data(i3 + 1) = data(i3rev + 1)

                         data(i3rev) = tempr
                         data(i3rev + 1) = tempi

                      enddo

                   enddo

                endif

                ibit = ip2 / 2

     1          continue
                if ( ( ibit .ge. ip1 ) .and. ( i2rev .gt. ibit ) ) then

                   i2rev = i2rev - ibit
                   ibit = ibit / 2
                   go to 1

                endif
                i2rev = i2rev + ibit

             enddo

             ifp1 = ip1

     2       continue
             if ( ifp1 .lt. ip2 ) then

                ifp2 = 2 * ifp1

                theta = isign * 3.14159265358979323846d0 * 2.0 /
         |           ( ifp2 / ip1 )
                wpr = -2.0 * sin( 0.5d0 * theta )**2
                wpi = sin( theta )

                wr = 1.0
                wi = 0.0

                do i3 = 1, ifp1, ip1

                   do i1 = i3, i3 + ip1 - 2, 2

                      do i2 = i1, ip3, ifp2

                         k1 = i2
                         k2 = k1 + ifp1

                         tempr = ( wr * data(k2) ) - ( wi * data(k2 + 1) )
                         tempi = ( wr * data(k2 + 1) ) +
         |                    ( wi * data(k2) )

                         data(k2) = data(k1) - tempr
                         data(k2 + 1) = data(k1 + 1) - tempi

                         data(k1) = data(k1) + tempr
                         data(k1 + 1) = data(k1 + 1) + tempi

                      enddo

                   enddo
                   wtemp = wr
                   wr = ( wr * wpr ) - ( wi * wpi ) + wr
                   wi = ( wi * wpr ) + ( wtemp * wpi ) + wi

                enddo

                ifp1 = ifp2

                go to 2

             endif

             nprev = n * nprev

          enddo

          return

          end
    следующий фpагмент (4)|пpедыдущий фpагмент (2)

    - Various nice sources (2:5030/84) ----------------------------- NICE.SOURCES -
     Msg  : 3442 of 3451
     From : Serguei Shtyliov                    2:5020/157.59   03 Jul 97  08:44:07
     To   : Stas Wilf                                           04 Jul 97  19:58:19
     Subj : FFT and SB
    -------------------------------------------------------------------------------
    Hail thou o Stas!

    Once upon a time thou didst write unto All:

     SW>   Hет ли у кого:

     SW>     Исходников FFT на Pascal

       Вот, вроде (не проверял):

    - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - -
    {$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.
    - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - - 8<- - - - -

    следующий фpагмент (5)|пpедыдущий фpагмент (3)
    - Various nice sources (2:5030/84) ----------------------------- NICE.SOURCES -
     Msg  : 8092 of 8108
     From : Vladimir L. Vasilevskij             2:5020/279.31   22 Nov 97  14:23:12
     To   : Andrew Kascheev                                     01 Dec 97  20:19:56
     Subj : FFT
    -------------------------------------------------------------------------------
    Hello Andrew!

     AK> У кого-ньть есть пpога, pеализующая FFT (Быстpое Пpеобpазование Фуpье)
     AK> на Pascal/C?

     В очередной раз....

     БПФ с замещением имени Кули-Тьюки по основанию 2.

    #include <math.h>

    #define  PI         3.1415926535
    #define  FFT        -1
    #define  REV_FFT    1

    void fft(float * x,float * y,int m,int d)   //  БПФ  по основанию 2

    // m - порядок (2**m точек)
    // d=1 - OБПФ, d=-1 - БПФ

    {

     int n,l,e,f,i,j,o,o1,j1,i1,k;

     float u,v,z,c,s,p,q,r,t,w,a;

     n=1<<m;

     for(l=1;l<=m;l++)
     {

     u=1;
     v=0;
     e=1<<(m-l+1);
     f=e/2;

     z=PI/f;

     c=cos(z);
     s=sin(z);

     if(d==FFT) s=-s;

     for(j=1;j<=f;j++)
      {
      for(i=j;i<=n;i+=e)
       {
       o=i+f-1;
       o1=i-1;
       p=x[o1]+x[o];
       r=x[o1]-x[o];
       q=y[o1]+y[o];
       t=y[o1]-y[o];
       x[o]=r*u-t*v;
       y[o]=t*u+r*v;
       x[o1]=p;
       y[o1]=q;
       }
      w=u*c-v*s;
      v=v*c+u*s;
      u=w;
      }
     }

     j=1;

     for(i=1;i<n;i++)
      {
      if(i<j)
       {
       j1=j-1;
       i1=i-1;
       p=x[j1];
       q=y[j1];
       x[j1]=x[i1];
       y[j1]=y[i1];
       x[i1]=p;
       y[i1]=q;
       }
      k=n/2;
      while(k<j)
       {
       j=j-k;
       k=k/2;
       }
      j+=k;
      }

     if(d==FFT) return;
      a=1.0/n;
      for(k=0;k<n;k++)
       {
       x[k]*=a;
       y[k]*=a;
       }
      return;
    }
    следующий фpагмент (6)|пpедыдущий фpагмент (4)
    - [101] Digital sound and music problems (2:5030/84) --------------- RU.STRACK -
     Msg  : 12 of 14
     From : Dmitry Niqiforoff                   2:5057/3        29 Jan 94  13:22:00
     To   : All
     Subj : Spectrum Analizer for GUS...
    --------------------------------------------------------------------------------
    .MSGID: 2:5057/3@fidonet 2d4a63d0
    * Crossposted in RU.STRACK
    * Crossposted in FIDO7.SOUND
    Hello All!

      Вот посмотрел тут я анализатор частот для GUSа, который в реальном времени
    рисует на экране столбики по частотам. Сколько отслеживается полос я даже
    сказать не берусь - что-то около 100, на частоте 44100 Hz. В докции есть
    кой-какое описание того, как это делается, и мне кажется многим это будет
    полезно, поэтому помещаю фрагмент этой документации здесь.

    ------------------------------------------------------------------------------
    ***************
    The FFT:

    The FFT is based on the Discrete Fourier Transform which is a discrete
    extension of the Fourier integral, an integral which can be interpreted as
    the continuous representation of the spectral components of any given
    function.  The Fourier integral itself derives from the Fourier series which
    is capable of describing any function over an interval in terms of a
    *unique* infinite summation of sines and cosines whose periods are integer
    multiples of each other.  Well, that was a mouthful but basically it is
    possible to describe a function, and sampled data over a range constitutes
    a function, as a summation of sine waves.  The mathematics of all this was
    worked out by the mega-geniuses of the past(Fourier,Laplace, and Euler)
    and boils down to this:
            1  N-1       jk
    X(j) =  -  SUM  x(k)W
            N  k=0

    where,
    N = number of points

    W = exp(2*pi*i)/N

    X(j) - N complex transformed points

    x(k) - N complex samples

    The points which are graphed on the screen are merely |X(j)| which come
    from the sample points x(k).  Since we are dealing with real data the
    complex parts of the samples are just set to zero.

    Now, as you can see the number of operations that must be gone through
    is proportional to N^2, N summations by N multiply-adds.  Up until 1963
    this is the way spectrums were done until Cooley-Tukey came up with a faster
    alternative that better utilizes intermediate results.  Their original
    paper described a method which worked with any number of samples so long
    as N was not prime and came to be know as the Fast Fourier Transform.  A
    special case arises when the number of samples is a power of 2(or a power of
    5 or any number for that matter). Spect uses just such an algorithm
    and runs in a time proportional to N*log<base 2>(N).  The time savings over
    N^2 is enormous.

    Its difficult to describe in words how the FFT algorithm utilizes
    intermediate results to reduce computation time.  Basically, it amounts
    to splitting up the set of samples into two odd and even sets, applying
    the summation detailed above except that exponents are all multiplied by two.
    With the resulting two sets any desired X(j) can be computed by adding an
    element from the even results with an appropriate element from the set of odd
    results multiplied by W raised to a fractional power depending on which X(j)
    is being computed.  The neat part is that this procedure can be applied
    recursively up to log2(N) times if N is a power of 2 until we produce N sets
    of one element each which are nothing more than the sample data themselves.
    Then the recursion stops and we can work our way backwards using the simple
    combining rules to produce an X(j) using only big O log2(n) operations to
    do so.  Over N X(j), this is proportional to N*log2(n) total compute time.

    The algorithm used by Spect I owe to who knows who but it dispenses with
    the stack pushing overhead involved in recursion and utilizes the patterns
    evident in doing all the substitutions of the combining rules into themselves
    that was affected by the recursion anyway.  The only draw back is that the
    results are in a scrambled order but a simple table lookup alleviates the
    headache of unscrambling them.  The code was written using 32bit integer
    instructions, utilizes table lookups wherever possible, and squeezes every
    last CPU cycle out of the i386.
    следующий фpагмент (7)|пpедыдущий фpагмент (5)
    - [42] Various nice sources (2:5030/84) ------------------------- NICE.SOURCES -
     Msg  : 39 of 44
     From : Serguey Zefirov                     2:5020/509.601  30 Nov 95  10:02:00
     To   : Pavel Shpin
     Subj : Анализатоp спектpа need !
    --------------------------------------------------------------------------------
    Zdorovenki bulji,(Hi! in other words) Pavel!
     Wednesday November 29 1995 17:07. Тогда-то Pavel Shpin и написал All буквально
    следующее:

     PS>     Пpиспичило вот мне - сигнал подается на вход SB.
     PS>  Далее надо по пpинципу пpогpаммы VPREG показывать на экpане
     PS>  каpтинку - pаскладку по спектpам. В общемм нужен алгоpитм
     PS>  Subj'а ....

    Беpете, и делаете FFT - Fast Fourier Transform.
    Вот тебе:

    -------------------------------8<------------------------------------
    #include <io.h>
    #include <mem.h>
    #include <math.h>
    #include <conio.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include <graphics.h>

    #define DISCR   11

    #define REDUCEDFOURIER 1
    #define REDUCELOW 0

    typedef unsigned short ushort;
    typedef unsigned char  uchar;

    struct complex
    {
        double re;
        double im;
    };

    const double PI = 3.14159265358979323846;

    complex *x;
    int n2 = DISCR;
    int nd = 1 << DISCR;
    int np = 1 << (DISCR - 1);

    inline double cabs( complex &num )
    {
        return sqrt( num.re*num.re + num.im*num.im );
    }

    // a - массив комплексных отсчетов
    // r - степень 2-ки. 1 << r - количество отсчетов
    // dir - напpавление пpеобpазования. 1-пpямое (спектp), -1-обpатное(тоже //
    спектp ;)
    // После обpаботки чеpез easy надо сделать scale, чтобы ноpмализовать амплитуды
    // Для частоты f и частоты семплиpования Fs амплитуда будет с индексом
    // (int)f*(1 << r)/Fs. Пpи этом надо учитывать, что полезными будут только //
    частоты вплоть до Fs/2. Выше - их зеpкальное отpажение. Вpоде все.
    // Если что - адpес знаете. ;)
    void easy( complex *a, int r , int dir)
    {
        ushort  n, lim1, lim2, i, j, l, m, k;
        double  arg, cs, si, cstep, sstep, cs1, si1;
        complex A, B, C;

        n = 1 << r;         // n - число отсчетов
        lim1 = n - 1;
        lim2 = n/2;

    // Inplace shuffle of data begins

        for( j = 1, i = 1; i <= lim1; i++ )
            {
            if( i < j )
            {
                A      = a[j-1];
                a[j-1] = a[i-1];
                a[i-1] = A;
            }

            l = lim2;
            while( l < j )
            {
                j -= l;
                l /= 2;
            }

            j += l;
            }

    // Shuffle complete
    // Transform by decimation in time begins

        for( i = 0; i < r; i++ )
            {
            lim1 = 1 << i;
            lim2 = 1 << (r-i-1);
            arg  = 2 * PI * lim2/n;
            cs = 1.0;
            si = 0.0;
     cstep = cos( dir*arg );
            sstep = sin( dir*arg );

            for( l = 0; l < lim2; l++ )
            {
                for( m = 0; m < lim1; m++ )
                {
                    k = m + l * 2 * lim1;
                    B = a[k];
                    C = a[k+lim1];
                    A.re =  C.re * cs + C.im * si;
                    A.im = -C.re * si + C.im * cs;
                    a[k].re = B.re + A.re;
                    a[k].im = B.im + A.im;
                    a[k+lim1].re = B.re - A.re;
                    a[k+lim1].im = B.im - A.im;
                    cs1 = cs*cstep - si*sstep;
                    si1 = si*cstep + cs*sstep;
                    cs = cs1;
                    si = si1;
                }
                cs = 1;
                si = 0;
            }
            }
    }

    void scale( complex *x, int r )
    {
        ushort n, i;

        n = 1 << r;

        for( i = 0; i < n; i++ )
        {
            x[i].re /= n;
            x[i].im /= n;
        }
    }

    -------------------------------8<------------------------------------
    следующий фpагмент (8)|пpедыдущий фpагмент (6)
    - [42] Various nice sources (2:5030/84) ------------------------- NICE.SOURCES -
     Msg  : 11 of 13
     From : Alexey Monastyrenko                 2:5030/303.8    23 Mar 96  21:08:00
     To   : All'ам
     Subj : БПФ на asm
    --------------------------------------------------------------------------------
     Приветствую , All'ам! ....

       Вот, давно обещаное БПФ на асме (большей частью) .

    === Cut === four_asm.pas

    {$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.
    следующий фpагмент (9)|пpедыдущий фpагмент (7)

    #include <io.h>
    #include <mem.h>
    #include <math.h>
    #include <conio.h>
    #include <stdio.h>
    #include <graphics.h>

    #define DISCR   11

    typedef unsigned short ushort;
    typedef unsigned char  uchar;

    struct complex
    {
        double re;
        double im;
    };

    const double PI = 3.14159265358979323846;

    void easy ( complex *a, int r );
    void scale( complex *x, int m );
    int  LoadData();
    void DrawSignal();
    void DrawFFT();

    complex *x;
    int n2 = DISCR;
    int nd = 1 << DISCR;
    int np = 1 << (DISCR - 1);

    inline double cabs( complex &num )
    {
        return sqrt( num.re*num.re + num.im*num.im );
    }

    void main( void )
    {
        if( LoadData() )
        {
            printf( "file error\n" );
            return;
        }

        int gd = DETECT, gm;
        initgraph( &gd, &gm, "" );

        DrawSignal();

        setcolor( 15 );
        settextjustify( RIGHT_TEXT, CENTER_TEXT );
        outtextxy( 639, 20, "FFT in progress" );

        easy( x, n2 );       // FFT
        scale( x, n2 );

        outtextxy( 639, 40, "Done!" );
        getch();

        DrawFFT();

        delete x;
        closegraph();
    }

    // a - массив комплексных отсчетов
    // r - степень 2-ки
    void easy( complex *a, int r )
    {
        ushort  n, lim1, lim2, i, j, l, m, k;
        double  arg, cs, si, cstep, sstep, cs1, si1;
        complex A, B, C;

        n = 1 << r;         // n - число отсчетов
        lim1 = n - 1;
        lim2 = n/2;

    // Inplace shuffle of data begins

        for( j = 1, i = 1; i <= lim1; i++ )
            {
            if( i < j )
            {
                A      = a[j-1];
                a[j-1] = a[i-1];
                a[i-1] = A;
            }

            l = lim2;
            while( l < j )
            {
                j -= l;
                l /= 2;
            }

            j += l;
            }

    // Shuffle complete
    // Transform by decimation in time begins

        for( i = 0; i < r; i++ )
            {
            lim1 = 1 << i;
            lim2 = 1 << (r-i-1);
            arg  = 2 * PI * lim2/n;
            cs = 1.0;
            si = 0.0;
            cstep = cos( arg );
            sstep = sin( arg );

            for( l = 0; l < lim2; l++ )
            {
                for( m = 0; m < lim1; m++ )
                {
                    k = m + l * 2 * lim1;
                    B = a[k];
                    C = a[k+lim1];
                    A.re =  C.re * cs + C.im * si;
                    A.im = -C.re * si + C.im * cs;
                    a[k].re = B.re + A.re;
                    a[k].im = B.im + A.im;
                    a[k+lim1].re = B.re - A.re;
                    a[k+lim1].im = B.im - A.im;
                    cs1 = cs*cstep - si*sstep;
                    si1 = si*cstep + cs*sstep;
                    cs = cs1;
                    si = si1;
                }
                cs = 1;
                si = 0;
            }
            }
    }

    void scale( complex *x, int r )
    {
        ushort n, i;

        n = 1 << r;

        for( i = 0; i < n; i++ )
        {
            x[i].re /= n;
            x[i].im /= n;
        }
    }

    int LoadData()
    {
        int f = _open( "dash.mrs", 0x8001 );
        int n;

        if( f < 5 )     return 1;

        _read( f, &n, sizeof( int ) );
        lseek( f, sizeof( int ) * n, SEEK_CUR );

        x = new complex[nd];
        double *src = new double[nd];

        _read( f, src, nd * sizeof( double ) );
        for( register i = 0; i < nd; i++ )
        {
            x[i].re = src[i];
            x[i].im = 0;
        }

        delete src;
        _close( f );

        return 0;
    }

    void DrawSignal()
    {
        double scale;
        int cy;

        setcolor( 4 );
        line( 0, 478, 639, 478 );           // horizontal axis

        setcolor( 1 );
        moveto( 0, 478 );

        scale = x[0].re;                    // scale input signal
        for( int i = 1; i < 512; i++ )
            if( x[i].re > scale )   scale = x[i].re;
        scale = 476.0 / scale;

        for( i = 0; i < 512; i++ )            // draw signal
        {
            cy = 478 - x[i].re*scale;
            lineto( i, cy );
            lineto( i+1, cy );
        }
    }

    void DrawFFT()
    {
        int y, base=0;

        line( 520,   0, 530,   0 );
        line( 520,  79, 530,  79 );
        line( 520, 159, 530, 159 );
        line( 520, 239, 530, 239 );
        line( 520, 319, 530, 319 );
        line( 520, 399, 530, 399 );
        line( 520, 479, 530, 479 );

        settextjustify( LEFT_TEXT, CENTER_TEXT );
        outtextxy( 535,   5, "0.0001" );
        outtextxy( 535,  79, "0.00001" );
        outtextxy( 535, 159, "0.000001" );
        outtextxy( 535, 239, "0.0000001" );
        outtextxy( 535, 319, "0.00000001" );
        outtextxy( 535, 399, "0.000000001" );
        outtextxy( 535, 475, "0.0000000001" );

        char c;

        setfillstyle( SOLID_FILL, 0 );
        setcolor( 10 );
        do {
            bar( 0, 250, 512, 479 );
            for( register i = 0; i < 512; i++ )
            {
                y = 480 - ((log10( cabs( x[i+base] ) + 1e-50 ) + 10.0) * 80.0);
                if( y > 479 )     y = 479;
                line( i, 479, i, y );
            }

            c = getch();
            switch( c )
            {
                case ',': if( base != 0 )           base -= 512; break;
                case '.': if( base != 512 )         base += 512;
            }

        } while( c != 27 );
    }

    следующий фpагмент (10)|пpедыдущий фpагмент (8)
    - [44] Nice sources discussion (2:5030/84) -------------------- NICE.SOURCES.D -
     Msg  : 24 of 24
     From : Anton Kondratyev                    2:5030/98.20    30 Dec 94  16:55:00
     To   : Andrey Elbakian
     Subj : FFT надо. (Пpеобpазования Фуpье однако)
    --------------------------------------------------------------------------------
        Hello, Andrey!

     On 29 Dec 94  02:07:00, Andrey Elbakian said to All
     (Re: FFT надо. (Пpеобpазования Фуpье однако))

     AE> Итак надобно мне для успешной дальнейшей pаботы алгоpитм быстpых
     AE> пpеобpазований Фуpье.

      [ преобразовано ]

     AE> Да, чуть не забыл, всю эту дpянь надо считать в pеальном вpемени с
     AE> пеpекpывающимися эпохами анализа для того чтобы можно было отследить
     AE> изменения спектpа (пpоцесс не стационаpный). Тоесть считать надо очень
     AE> быстpо... А там еще минимум 16 каналов...
     AE> Если у кого есть любая инфоpмация, а еще лучше алгоpитмы, не дайте
     AE> погибнуть.
     AE> Заpанее спасибо.

      Итак . Программа вычисления коэффициентов дискретного преобразования Фурье.
      (на Фортране, правда)

       SUBROUTINE TRFUR (X,A,B,N,G)
       DIMENSION X(N),A(N),B(N)
       G=0
       A1=0
       AN=N
       DO1 I=1,N
       G=G+X(I)
    1  A1=A1+(-1.)**I*X(I)
       G=G/AN
       A(N/2)=A1/AN
       M=N/2-1
       DO2 I=1,M
       AI=I
       A(I)=0
       B(I)=0
       DO3 J=1,N
       AJ=J
       T=6.2832*AI*AJ/AN
       C=COS(T)
       S=SIN(T)
       A(I)=A(I)+X(J)*C
    3  B(I)=B(I)+X(J)*S
       A(I)=A(I)*2./AN
    2  B(I)=B(I)*2./AN
       RETURN
       END

      X - вещественный массив обрабатываемых чисел
      A,B - массивы вещественных и мнимых частей ДПФ
      N - длина входного массива
      G - постоянная составляющая ДПФ.
    Hашел я это в книге С.И.Баскаков "Радиотехнические цепи и сигналы",
           М,"Высшая школа",1988
    Это не БПФ, но затруднений не должно быть, берешь N=2^p и делишь там все
    в соответствии с формулой. А вообще интересная проблема!

    Всего 9 фpагмент(а/ов) |пpедыдущий фpагмент (9)

Если вы хотите дополнить FAQ - пожалуйста пишите. design/collection/some content by Frog, DEMO DESIGN FAQ (C) Realm Of Illusion 1994-2000, При перепечатке материалов этой страницы пожалуйста ссылайтесь на источник: "DEMO.DESIGN FAQ, http://www.enlight.ru/demo/faq".