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".