следующий ф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 и делишь там все
в соответствии с формулой. А вообще интересная проблема!