logo elektroda
logo elektroda
X
logo elektroda
REKLAMA
REKLAMA
Adblock/uBlockOrigin/AdGuard mogą powodować znikanie niektórych postów z powodu nowej reguły.

Algorytm FFT. Dla większej liczby próbek pojawia się błąd.

Prymulka 08 Paź 2006 16:59 8318 9
REKLAMA
  • #1 3097397
    Prymulka
    Poziom 18  
    Witam

    Wklepałem sobie program to Matlaba z książki Zielińskiego "CPS od teorii do zastosowań". Wygląda on tak:

    %GENERACJA SYGNAŁU 
    N=8; %liczba próbek sygnału 
    x=0:N-1; %przykładowe wartości próbek 
    typBitReverse=1; 
    typFFT=1; 
    
    %PRZESTAWIENIE KOLEJNOŚCI PRÓBEK: wersja 1-wolna 
    ...........
    
    %WŁAŚCIWE FFT - wersja 1 - wolniejsza 
    if (typFFT==1) 
      for e=1:log2(N) 
        SM = 2^(e-1); 
        LB = N/(2^e); 
        LMB = 2^(e-1); 
        OMB = 2^e; 
        W = exp(-j*2*pi/2^e); 
        for b = 1 : LB 
          for m = 1 : LMB 
            g = (b-1)*OMB + m; 
            d = g + SM; 
            xgora = x(g); 
            xdol = x(d)*W^(m-1); 
            x(g) = xgora + xdol; 
            x(d) = xgora - xdol; 
         end 
      end 
     end 
    end 


    No i niby wszystko działa jak należy. Problem się zaczyna gdy wstawię większą liczbępróbek np. N = 128. W linii
    xgora = x(g);
    wyskakuje mio bład że g jest za duze (tzn wieksze niz maksymalny x czyli 127). No i jak sobie wszystko w głowie popodliczam to rzeczywiscie tak jest. Tylko nie mam pojecia co zrobic zeby tego bledu uniknac :/.
    Czemu ten algorytm działa tylko dla N<=8?

    W załączniku znajduje sie caly plik z matlaba. W tym momencie mam ustawione N = 64 więc wyskoczy blad po odpaleniu pliku 'radix.m'
  • REKLAMA
  • #2 3103426
    Xitami
    Poziom 29  
    A czemu SM == LMB

    x[0..n-1] a kiedy "d" lub "g" będzie równe zero?
  • #3 3104402
    Prymulka
    Poziom 18  
    Witam

    Dzieki za zainteresowanie
    SM - szerokość motylka
    LMB - liczba motylkow w bloku

    "d" i "g" nie osiągną zero ponieważ z tego co się orientuję matlab indeksuje tablice od "1". Tak więc pierwsza wartość w tabeli x[1] = 0
  • REKLAMA
  • #4 3105575
    Xitami
    Poziom 29  
    W kodzie który zamieściłeś, zawsze SM==LMB==2^(e-1), no to jakby o jedną zmienną za dużo, albo cosik źle przepisał kolega.
    X=0:N-1 czyli sygnał znajduje się w x od indeksu 0 do indeksu N-1, piszesz że błąd pojawia się gdy „g”>127 (dla N=128) i tu jakby rozumiem, czyli minimalna wartość indeksu musi być wtedy równa zero.
    A na moje oko „g” biega od jedynki, no to może powinno być „x=1:N”? (jeżeli jest to deklaracja tablicy, bo nie znam MatLaba.)
  • #5 3106901
    Prymulka
    Poziom 18  
    %WŁAŚCIWE FFT - wersja 2 - szybsza
    if (typFFT==2)
        for e = 1: log2(N)
            L = 2^e;
            M = 2^(e-1);
            Wi = 1;
            W = cos(2*pi/L)-j*sin(2*pi/L);
            for m = 1 : M
                for g = m : L : N
                    d = g + M;
                    T = x(d)*Wi;
                    x(d) = x(g) - T;
                    x(g) = x(g) + T;
                end
                Wi=Wi*W;
            end
        end
    end


    Tamta zmienna jest celowo napisana dwa razy ponieważ jest to przykład mający dobrze obrazować działanie FFT. Ten powyżej jest juz szybszy. W linii
    T = x(d)*Wi;
    wyskakuje mi ten sam błąd. Zmiana x nic nie daje.
  • Pomocny post
    #6 3108633
    Xitami
    Poziom 29  
    Algorytm II jest poprawny!!!
    Porównałem go nawet z algorytmem ze str. 102 w V wydaniu „Metody numeryczne” Z. Fortuna, B. Macukowski, J. Wąsowski, w książce tej jest błąd, powinno być „W = CMPLX(COS(ANG), -SIN(ANG))” (brakowało znaku "-").
    Indeksy przebiegają wartości od jeden do N (nie do N-1) i tak właśnie należy zadeklarować i zainicjować tablicę x().
  • REKLAMA
  • #7 3126762
    Prymulka
    Poziom 18  
    Miałem w programie czeski bład :) Ale mniejsza o to. Wlasnie robie projekt na procku sygnałowym odnośnie analizy widma. LCD graficzny juz chodzi a z fft narazie napisalem motylka. Teraz chcę przerobić powyzszy algorytm aby obliczal analize tylko dla sygnałow rzeczywistych. Jak najprosciej sie za to zabrac?
  • #8 3127477
    Xitami
    Poziom 29  
    ...mniejsza o to...” no pięknie. ;-)

    Zobacz np. 12 rozdział w www.dspguide.com.

    Albo
    Próbki umieść w tablicy X w następującej kolejności: 0,2,4,...,n-2,1,3,...n-1 (parzyste w części RE [pierwsza połowa tablicy], a nieparzyste w części IM [druga połowa tablicy]).
    Porządek taki można uzyskać np. dzięki następującej procedurze (Pascal):
    procedure oddsort(ncap:integer; var y:longvector);
    //w miejscu sortuje ciąg, tab by uzyskać 0,2,4,...,n-2,1,3,...n-1
    var
      tcap, t, j, k:integer;
      store:float;
    begin
      tcap:= ncap div 2;
      for j:=1 to tcap-1 do begin
        k:=j;
        if j>1 then
          repeat
            k:=(k*2) mod (2*tcap-1);
          until k<=j;
        if k=j then begin
          store:=y[j];
          t:=j;
          repeat
            k:=t;
            t:=(2*k) mod (2*tcap-1);
            y[k]:=y[t]
          until t=j;
          y[k]:=store
        end
      end
    end;
    Ale zgrabniej będzie chyba już u źródła utykać je w takim porządku.

    Następująca procedura załatwi „kompaktowe” FFT:
    procedure CompactRealFFT(var x:longvector; ncap, g:integer);             // ncap=2^g
    var
      t, tcap:integer;
      xReven, xRodd, xIeven, xIodd, store:float;
      theta, increment, sine, cosine:float;
    begin
      tcap:= ncap div 2;
      increment:= pi/tcap;
      theta:= 0;
      oddsort(ncap, x);                         // jak pisałem wcześniej
      g:=g-1;
    
      base2fft(x, tcap, g);                     // o tym niżej
    
      for t:=1 to tcap div 2-1 do begin        
        theta:= theta + increment;             
        cosine:= cos(theta);                    
        sine:=   sin(theta);                    
        xReven:= (x[t     ] + x[tcap-t])*0.5;
        xRodd := (x[t     ] - x[tcap-t])*0.5;
        xIeven:= (x[t+tcap] + x[ncap-t])*0.5;
        xIodd := (x[t+tcap] - x[ncap-t])*0.5;
        store:=               cosine*xIeven - sine*xRodd;
        x[t     ]:=  xReven + store;
        x[tcap-t]:=  xReven - store;
        store:=               cosine*xRodd  + sine*xIeven;
        x[t+tcap]:=  xIodd  - store;
        x[ncap-t]:= -xIodd  - store;
      end;
    
      store:=x[0];
      x[0]:=x[0]+x[tcap];
      x[tcap]:=store - x[tcap]
    end;
    base2fft(x, tcap, g); pierwsza połowa wektora to część rzeczywista, druga urojona, tcap=2^g

    D.S.G. POLLOCK, (1999), A Handbook of Time-Series Analysis, Signal Processing and Dynamics, Queen Mary and Westeld College The University of London UK
  • REKLAMA
  • #9 3129461
    Prymulka
    Poziom 18  
    Dzieki za odpowiedz
    Posegregowalem próbki w tablicy na parzyste i nieparzyste. Teraz pytanie odnosnie tego base2fft. Jakie fft ona liczy. Bo przeciez ponizej jest fft. Przeanalizowalem przyklad w basicu z dspguide i tam jest najpierw policzone fft z czesci parzystej a pozniej reszta
  • #10 3133349
    Xitami
    Poziom 29  
    base2fft(x, tcap, g) to zwykłe FFT, z tym że w wektorze x[0..tcap*2-1], początek [0..tcap-1] to część rzeczywista, a koniec [tcap..2*tcap-1] część urojona, tcap==2^g. Dodam trochę kodu (z tego samego żródła):
     procedure base2fft(var y:longvector; tcap, g:integer);
    var
      a, c, t, i, j, k, l, P,Q : integer;
      W, theta, sine, cosine,
      yR, yI : float;
    begin
      w:=twopi/tcap;
      p:=tcap;
      c:=0;
      for l:=1 to g do begin
        a:=0; t:=0;
        if l=1 then q:=1 else q:=q*2;
        p:=p div 2;
        for i:=0 to q-1 do begin
          for k:=0 to p-1 do begin
            t:=a+c;
            yr:=y[t     ];
            yi:=y[t+tcap];
            y[t       ]:= yr + y[t+p     ];
            y[t+tcap  ]:= yi + y[t+p+tcap];
            y[t+p     ]:= yr - y[t+p     ];
            y[t+p+tcap]:= yi - y[t+p+tcap];
            if l<g then begin
              theta:=w*((c*q) mod tcap);
              cosine:= cos(theta);         // c*q indeks do ew. tablic
              sine:=   sin(theta);
              yr:= y[t+p     ];
              yi:= y[t+p+tcap];
              y[t+p     ]:= yr*cosine + yi*sine;
              y[t+p+tcap]:= yi*cosine - yr*sine;
            end;
            c:=c+1;
          end;
          c:=0;
          a:=a+2*p
        end
      end;
    
      for j:=1 to tcap-2 do begin
        t:=bitreverse(j,g);
        if t>j then begin
          yr:=y[t];
          yi:=y[t+tcap];
          y[t]:=y[j];
          y[t+tcap]:=y[j+tcap];
          y[j]:=yr;
          y[j+tcap]:=yi;
        end
      end
    end; {base2fft}
    
    function bitreverse(j,g:longint):longint;
    var t,tl,r,l:longint;
    begin
      t:=0; r:=j;
      for l:=1 to g do begin
        tl:=r mod 2;
        t:=t*2+tl;
        r:=r div 2
      end;
      result:=t
    end;
    W zasadzie zwykłe DIF FFT, choć sortowanie "bitreverse" zrobiłbym raczej "normalnie".
REKLAMA