wyrysowac wykres

Witam!! Jak wyrysowac wykres widma amplitudowego sygnału? Cała funkcja jest obliczana za pomocą FFT i do tego na ekranie musi mi jesszcze wyświetlić wykres. Poniżej cały schemat obliczenia, chociaz nie wiem czy do końca jest poprawny, bo może zawierać błedy, mógłby go ktoś przy okazji sprawdzić??

const
NN = 256; //liczba probek
LPI = 3.14159265358979323846;

procedure TForm1.fourier(Sender: TObject; X,Y: array of Double);
//procedure TForm1.fourier(Sender: TObject; X,Y: array of 

Double):[ewentualnie_typ_zwracany]; // jeli co zwraca??


  var
    a,z,licznik : Integer;
    b,c,d,e,g,L,M,mm,i: Integer;
    W_r,W_i,Wi_r,Wi_i,T_r,T_i: Double;
    //licz[1024]: real;
    T : array [] of real;

   licznik:= 0;
   a:=1;


//-----------------------------------------przestawianie elementow

   for b:=1  to b<NN  do

      begin
        if b<a  then

                T:=X[a-1];
                X[a-1]:=X[b-1];
                X[b-1]:=T;
            end;
       
        c:=NN/2;

        while c<a do
              begin
                 a:=a-c;
                 c:=c/2;
              end;
       
        a:=a+c;  //??

      b:=b+1;
      end;

//--------------------------------------FFT

    for e:=0 to e<NN do
    begin
         Y[e]:=0;
    
    e:=e+1;
    end


    for e:=1 to e<=9 do
       begin
           L:= 1<<e;  // z tego co obczaiłem pow(1,e) niewiem trzeba by to 

jako napisac procedurke 
                     //albo jest może gotowa implementacja w delphi
           M := L/2;
           Wi_r:=1;
           Wi_i:=0;

           W_r= cos(LPI/M);
           W_i=-sin(LPI/M);
      
           for mm:=1 to mm<M+1 do		// jak w pocie można 

chyba jako zlikwidować M+1, NN+1
              begin

                  for g:=mm to g<NN+1 do 
                      begin
                           d:= g+M;
                           
                           //T=x(d)*Wi
                            
                           T_r := X[d-1]*Wi_r-Y[d-1]*Wi_i;
                           T_i := X[d-1]*Wi_i+Y[d-1]*Wi_r;

                            // x(d)=x(g)-T

                           X[d-1] := X[g-1]-T_r;
                           Y[d-1] := Y[g-1]-T_i;
                          
                           //x(g)=x(g)+T
                           
                           X[g-1] := X[g-1]+T_r;
                           Y[g-1] := Y[g-1]+T_i;
		     
                      g:=g+L;   
                      end;

              Wi_r:=Wi_r*W_r-Wi_i*W_i;
              Wi_i:=Wi_r*W_i+Wi_i*W_r;
	      
	      mm:=mm+1;
              end;
    e:=e+1;
    end;


end;

tak na pierwszy rzut oka:

              Wi_r:=Wi_r*W_r-Wi_i*W_i;
              Wi_i:=Wi_r*W_i+Wi_i*W_r;

pamiętasz?

niebardzo, o co z tym chodzi? Aha, jeszcze jedno pytanko, jak wygenerować sygnał rzeczywisty w delphi, aby został poddany analizie.

Myślałem, że już raz to pzerabiałeś
Wi_r:=Wi_rW_r-Wi_iW_i;
Wi_i:=Wi_rW_i+Wi_iW_r;
A widzisz różnicę z:
T_r:= Wi_rW_r-Wi_iW_i;
Wi_i:=Wi_rW_i+Wi_iW_r;
Wi_r:=T_r;

Rozumiem, że X to część rzeczywista, a Y to urojona. Pamiętaj o wyzerowaniu Y przed liczeniem FFT. 
Czemu T jest tablicą? 
Wykres to sqrt(sqr(X[i]) + sqr(Y[i]) dla i:=0..N/2

Dam Ci kawałek sprawdzonego kodu, indeksy liczone od zera do N-1. M=log2(N), albo inaczej 2^N=M

lm:=0;
j:=nd2;
for i :=1 to n-2 do begin
if i<j then begin
tr:=rex[j]; rex[j]:=rex[i]; rex[i]:=tr;
{ti:=imx[j]; imx[j]:=imx[i]; imx[i]:=ti;}
{^ rozkomentować dla sygnału zespolonego}
end;
k:=nd2;
while k<=j do begin
j:=j-k; k:=k div 2
end;
j:=j+k;
end;

for L:=1 to m do begin
le:=1 shl L; le2:=le div 2;
ur:=1; ui:=0;
sr:= cos(pi/le2);
si:=-sin(pi/le2);
for j:=1 to le2 do begin
i:=j-1;
while i<=nm1 do begin
lm:=lm+1;
ip:=i+le2;
tr:=rex[ip]ur - imx[ip]ui;
ti:=rex[ip]ui + imx[ip]ur;
rex[ip]:=rex[i] - tr;
imx[ip]:=imx[i] - ti;
rex[i] :=rex[i] + tr;
imx[i] :=imx[i] + ti;
i:=i+le
end;
tr:=ur;
ur:=trsr - uisi;
ui:=trsi + uisr;
end
end;

Myślałem, że już raz to pzerabiałeś
Wi_r:=Wi_rW_r-Wi_iW_i;
Wi_i:=Wi_rW_i+Wi_iW_r;
A widzisz różnicę z:
T_r:= Wi_rW_r-Wi_iW_i;
Wi_i:=Wi_rW_i+Wi_iW_r;
Wi_r:=T_r;

Rozumiem, że X to część rzeczywista, a Y to urojona. Pamiętaj o wyzerowaniu Y przed liczeniem FFT. 
Czemu T jest tablicą? 
Wykres to sqrt(sqr(X[i]) + sqr(Y[i]) dla i:=0..N/2

Dam Ci kawałek sprawdzonego kodu, indeksy liczone od zera do N-1. M=log2(N), albo inaczej 2^N=M

lm:=0;
j:=nd2;
for i :=1 to n-2 do begin
if i<j then begin
tr:=rex[j]; rex[j]:=rex[i]; rex[i]:=tr;
{ti:=imx[j]; imx[j]:=imx[i]; imx[i]:=ti;}
{^ rozkomentować dla sygnału zespolonego}
end;
k:=nd2;
while k<=j do begin
j:=j-k; k:=k div 2
end;
j:=j+k;
end;

for L:=1 to m do begin
le:=1 shl L; le2:=le div 2;
ur:=1; ui:=0;
sr:= cos(pi/le2);
si:=-sin(pi/le2);
for j:=1 to le2 do begin
i:=j-1;
while i<=nm1 do begin
lm:=lm+1;
ip:=i+le2;
tr:=rex[ip]ur - imx[ip]ui;
ti:=rex[ip]ui + imx[ip]ur;
rex[ip]:=rex[i] - tr;
imx[ip]:=imx[i] - ti;
rex[i] :=rex[i] + tr;
imx[i] :=imx[i] + ti;
i:=i+le
end;
tr:=ur;
ur:=trsr - uisi;
ui:=trsi + uisr;
end
end;

Pamiętam o wyzerowaniu. A tą część co mi podałeś to powinienem umieścić zamiast mojej części kodu? Jeszcze o coś chcĘ zapytać, umiałbyś napisać taki programik do obliczania widma amplitudowego dla sygnału rzeczywistego? Daj znać, nawet na maila: [email protected]

Przez pomyłkę poprzedni post wysłał mi się dwukrotnie.
Tak, to jest to co u Ciebie, zmienne inaczej się nazywają i tyle.

for i:=0 to n-1 do Y[i]:=0;
fft(X,Y);
for i:=0 to n div 2 do
    widmo[i]:=sqrt(sqr(X[i]) + sqr(Y[i])

(mam nadzieję że to co posłałem pocztą też się przydało)

właśnie to obczajam:) ten kod powinienem wstawić zaraz po kodzie wczytywania danych, bo mi wywala błędy jak tworze nową procedure?