Elektroda.pl
Elektroda.pl
X

Search our partners

Find the latest content on electronic components. Datasheets.com
Elektroda.pl
Please add exception to AdBlock for elektroda.pl.
If you watch the ads, you support portal and users.

Potrzebna pomoc w Dyskretnej Transformacie Fouriera

slawek55 25 Sep 2006 06:54 16330 45
  • #31
    slawek55
    Level 23  
    CZesc
    To ile tablic jest mi potrzebnych. Powiedzmy że mamay tablice spectrum po obliczeniach określoną spectrum[N] i jaki jeszcze będą potrzebne?

    Dodano po 53 [minuty]:

    Ale z tymi logarytmami to jest ne do mojego przypadku, czy tak. JEzeli mam zmienna spectrum typu uint16_t a ja zniej mam zrobić logarytm o wartości max 64. To co z tym.


    Jeszcze chciałbym poprosić o szrszy opis opadania słupków widma (jesli moge). Chodzi to ze w przewaniu nastepuje uaktualnienie opadania a w pętli głównej z tablicy spectrum robisz aktualną watość prążka. i czy do każdego prążka mam inna warośc czasu?
  • #32
    Xitami
    Level 29  
    Liczyłem logarytmy nie ze spektrum, a z wartości zespolonej (zauważ brak pierwiastkowania).
  • #33
    slawek55
    Level 23  
    Ok. rozumiem.
    tyle ze spectrum ma zakres od 0 do około 35000 a wynik ma być w zakresie 0 do 64. Czy to tak jest
  • #34
    slawek55
    Level 23  
    Mam prośbę nie wiem jak sobie poradzić z tym logarytmowaniem czy może mi ktoś pomóc napisać to.

    Druga rzecz która nie daje mi spokoju to obsługa słupków i pick-hold. Czy dla każdego słupka musi być osobna zmienna przechowująca czas opóźnienia, troche brakuje mi pamięci.
  • #35
    Prymulka
    Level 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
    if (typBitReverse==1)
    MSB=log2(8);
    for n=0:N-1;
    ncopy=n;
    nr=0;
    for m=1:MSB
    if(rem(n,2)==0)
    n=n/2;
    else
    nr=nr+2^(MSB-m);
    n=(n-1)/2;
    end
    end
    y(nr+1)=x(ncopy+1);
    end
    x=y;
    end

    %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?
  • #36
    Xitami
    Level 29  
    Cześć, nie odpowiadałem ostatnio bom miał dość pijany koniec tygodnia. A to kolegę żona odeszła, a to prokuratorowi coś się przypomniało (trzeba było uczcić że niepotrzebnie), a to znowu mnie dopadły koszmary ;-).
    Na otrzeźwienie poczytałem sobie „Cyberiadę”. Nie do końca podziałała, ale może dla tego że, w lodówce było jeszcze trochę płynnego jęczmienia?
    Szkoda, że Lem’cio ostatnio uciekł w filosofię (i wspominki) z której tak pięknie szydził. Choć może w ten sposób zamknął jakoś swoje dzieło, by mniej nam było Go brak.

    Ale wróćmy do dziesiątych pana Bell’a.
    Maksimum w twoim widmie to X==65’614, trochę mnie to martwi, bo to więcej niż Atmelkowe MAXINT (czyli 2^(16-1) ale jest to może wynikiem sumowania oktawowego, no cóż, nasze ucho nie zna się na matematyce, ale zna się (ucho) na wynikach doświadczeń opisanych przez kolegę Darwina.

    Zakładam, że prążek jest 16 bitową liczbą bez znaku („unsigned” co troszkę komplikuje zabawę).
    I ja i Shag’u kombinowaliśmy z logarytmowaniem przy podstawie 2, proponowaliśmy podnoszenie do 4 potęgi albo do 2 (czyli brak pierwiastka z modułu, o czym pisałem) ale przecie to nie dodaje nowej informacji. A zróbmy sobie tabelkę (LUT z angielska) np. taką:
    Code:
    flash const unsigned int t[]=(
    
    33622,29966,26707,23803,21214,18907,16851,15019,13386,11930,
    10633, 9476, 8446, 7527, 6709, 5979, 5329, 4750, 4233, 3773,
     3363, 2997, 2671, 2381, 2122, 1891, 1686, 1502, 1339, 1193,
     1064,  948,  845,  753,  671,  598,  533,  475,  424,  378,
      337,  300,  268,  239,  213,  190,  169,  151,  134,  120,
      107,   95,   85,   76,   68,   60,   54,   48,   43,   38,
       34,   30,   27,   24,   24);

    wtedy dla zadanego S gdy zawołamy:
    Code:
    j=0;
    
    if (S<t[64]) j=65;
    else while (t[j]>S)
          j += 1;
    Wartość „1-j” odpowiada zaokrąglonej wartości 20*log(S/X) (X==35’614)
    Tabelkę „t[]” stworzyłem tak:
    Code:
      p=69;      //poprzednia wartość
    
      for (i=X+10000; i>=23; j--) {                 // 23 bo 64 ;-)
        db=round(20*log10(i/X));
        if (db != p) {
          t[-db]=i+1;       //nie pamiętam czemu napisałem „+1”
          p=db;
        }
      }
    Chyżej można przeglądać tabelkę „bisekcyjnie”, może warto?

    Dygresja:
    Code:
    w C pisze się i kompiluje się efektywniej, lecz Pascal daje „magiczną” możliwość nadzorowania (nie lubię w sumie angielskiego słowa „controla”) zakresu zmiennych, indeksów itp, zauważyłem wiele „postuf” na forum (tym czy innym) wynikających tylko z tego, że np. indeks zagalopował[i] się zbyt daleko (czy blisko), albo że zmienna po dodaniu dodatniej wartości staje się ujemna. C daje nam <u>asercje</u> Fajna mechanika, lecz polega na „świadomości” programisty. W innym języku też to można zrobić, lecz C potrafi jakoś gładziej wygonić warunki z kompilowanego (na gotowo) kodu (w innych językach trzeba stosować warunkową kompilację). Polecam młodszym kolegom zaznajomienie się ze z tym sposobem. 
    
    Pytanko: jak zgrabniej przetłumaczyć słowo „assert”? „Asercja” jakoś drapie mnie w ucho, a litera „ł” nie może wystąpić w identyfkatorze.
    Druga dygresja:
    Code:
    Język C, to już w zasadzie synonim zwięzłości (patrz np.:[url]http://www.ioccc.org/[/url]) a bywa że przegrywa pod tym względem ze staruszkiem FORTRANem.
  • #37
    shg
    Level 35  
    Powiem tak: diablo sprytnie, czyli szybko i dokładnie :]
    Dla 64 punktów bedzie jak znalazł. Powinno działać szybciej niż operacje na 64 bitach.
    Nie wiem, co na to kompilator, ale samą pętlę można zrealizować w 4 instrukcjach + skok:
    Code:
    pyntla:
    
    LD Rx, Z+ ; zaladowanie dwoch kolejnych bajtow z tabeli
    LD Ry, Z+
    CP Rz, Rx ; porownanie z amplituda sygnalu
    CPC Rw, Ry
    BRCS pyntla

    Rx, Ry - do tych rejestrow ładowana jest wartość z tabeli.
    Rz, Rw - w tych rejestrach przechowywana jest aktualna amplituda sygnału
    x i z - bajt mniej znaczący, y i w - bajt bardziej znaczący.
    Do rejestru Z przed rozpoczęciem wykonywania pętli ładowany jest wskaźnik na początek tablicy (w RAM), w sumie można by ją równie dobrze umieścić w pamięci programu, wtedy trzeba użyć LPM zamiast LD. odczyt z RAM jest szybszy, ale tablica zajmuje cenną pamięć.
    Z warunkiem BRCS tak sobie strzeliłem tylko, może być inny (nigdy nie lubiałem skoków warunkowych :P), w każdym razie chodzi o skok gdy Rx:Ry będzie większe od Rz:Rw.
    Oczywiście należy dopilnować żeby wartości w tablicy miały odpowiedni porządek bajtów (little-/big-endian) w razie potrzeby zamienić miejscami dwie pierwsze instrukcje.
    Potem od tego co było w Z odejmujemy adres początku tablicy i to w zasadzie tyle.
    Żeby mieć odrazu wysokość słupka można odwrócić wyszukiwanie (i kolejność elementów w tablicy) w ten sposób z wartości Z otrzymamy odrazu wysokość słupka.

    A czy bisekcja będzie szybsza to nie wiem, przy większej tabeli na pewno, ale tu może okazać się, że overhead okaże się na tyle duży, że lepiej będzie użyć 'naiwnego' algorytmu przeszukiwania, czyli kolejno porównywać.

    Dla każdego słupka wartość opóźnienia jest stała, ale każdy musi mieć swoją zmienną (a najlepiej w jednej tablicy to trzymać) w której jest licznik czasu opóźnienia po którym markery wartości szczytowych zaczynają opadać.
    u mnie jest tak:
    volatile uint16_t bars[N_BANDS]; /* wartości poszczególnych słupków w skali liniowej */
    volatile uint8_t bands[N_BANDS]; /* wysokość słupków */
    volatile uint8_t peaks[N_BANDS]; /* pozycja wartości szczytowej */
    volatile uint8_t peak_delay[N_BANDS]; /* wyżej wymienione liczniki */
    volatile uint8_t leds[N_BANDS]; /* to co będzie wyświetlane, każdy bit to jedna dioda, taka bitmapa w sumie */

    Każdy słupek musi mieć osobny licznik czasu, bo inaczej zaczęły by opadać wszystkie w tym samym momencie i nie wiadomo w jakim, bo jakoś i tak by trzeba to zliczać, a jak jest jeden licznik to nie za bardzo się da (a właściwie to się da, tylko nie wiadomo, kiedy należało by go resetować).
    Opadanie samych słupków nie wymaga osobnych zmiennych, potrzebna jest tylko tablica, która przechowuje ich wysokość i operuje się bezpośrednio na tej tablicy zmniejszając co pewnien ustalony czas każdy z jej elementów o 1 (o ile jest co zmniejszać, czyli tylko gdy dany element jest większy od 0).
    To samo robimy z licznikami czasu dla każdego słupka - co jakiś czas zmniejszamy o 1 liczniki większe od zera.
    I to samo z wartościami szczytowymi, ale dodatkowy warunek - zmniejszamy je tylko gdy odpowiadające im liczniki czasu mają wartość równą zero.
    Liczniki czasu resetujemy (zapisujemy do nich żądaną wartość opóżnienia, nie zero) gdy wyliczona wysokość słupka jest większa od aktualnej pozycji markera wartości szczytowej.
    Podobnie ze słupkami, jeżeli ma działać powolne opadanie, to elementy tablicy przechowującej wysokości słupków aktualizujemy wyłącznie jeżeli wyliczona wysokość jest większa od tej w tablicy, jak nie jest, to ją ignorujemy.
  • #38
    slawek55
    Level 23  
    cześć.
    Czy mogę mieć do Ciebie prośbę?
    Nie mam takiego doświadczenia w pisaniu prog. i powiem szczerze że teraz już nic z tego nie rozumiem. Co mam zrobić z tym logarytmem. jak teraz to ubrać w funkcję (najlepiej), która zwróci mi wynik właściwy z tablicy spectrum.

    Chciałbym wrócic jeszcze do słupków i picków. Jak zamiast opadania pika srawić aby znikał. Codzi mi o to, ze widziałem w jednym sprzecie taki wskażnik, który też miał pick-hold, ale ten pik nie opadał tylko po pewnym czasie znikał i pojawiał sie nowy (po prostu pozbycie sie tylko opadania)

    Bardzo prosze o pomoc, wiem że mozesz powiedzieć skoro nie do końca wiesz o co tu chodzi to sie nie zabieraj za to, ale z drugiej strony chciałbym bardzo dokończyć sobie to i zobaczyć efekt ciężkiej pracy oraz przy okazji sie duzo nauczyć (z waszą pomocą).

    Sławek.
  • #39
    Xitami
    Level 29  
    Code:
     int decybele(unsigned int S){
    
        unsigned char j=0;
        else while (t[j]>S)
            j+=1;
        return(1-j);
    }
  • #40
    slawek55
    Level 23  
    Dzieki Xitami że odpowiedziałeś.



    A możesz mi powiedziec cos o tej tablity t[] i jej powstaniu?

    1. Skad liczba 69 oraz w petli for x+10000


    2. czy ta funkca to jakby odpowiednik wzoru 20*log10(x/35614)
    Poniewaz do logarytmu podstawiasz liczbe mniejsza od 1 to wyjdzxie ujemny a my operujemy na liczbach unsigned.
  • #41
    shg
    Level 35  
    Ze znikaniem tak samo jak z opadaniem - też trzeba liczniki mieć i to w takiej samej ilości. Różnice:
    Nie zmniejszamy co określony czas elementów tablicy przechowujących wartości szczytowe. Zmniejszamy tylko liczniki czasu.
    Podczas wyświetlania rysujemy tylko te punkty, których liczniki czasu są większe od zera.
  • #42
    Xitami
    Level 29  
    Przejrzałem i uładziłem.
    Code:
    flash const unsigned int t[64]={  // 0 .. –64 dB, tym sposobem można do –77dB
    
    33622,29966,26707,23803,21214,18907,16851,15019,
    13386,11930,10633, 9476, 8446, 7527, 6709, 5979,
     5329, 4750, 4233, 3773, 3363, 2997, 2671, 2381,
     2122, 1891, 1686, 1502, 1339, 1193, 1064,  948,
      845,  753,  671,  598,  533,  475,  424,  378,
      337,  300,  268,  239,  213,  190,  169,  151,
      134,  120,  107,   95,   85,   76,   68,   60,
       54,   48,   43,   38,   34,   30,   27,   24};

    signed char decybel(unsigned int s) {   // 20*log10(s/35614)
       char j=0;
       if (s < t[63]) j=64;
       else while (s < t[ j])
          j++;
       return(-j);    //jeżeli przeszkadzają ujemne zwróć 64-j
    }
    tabelka powstała tak:
    Code:
      p:=0;
    
      for i:=x downto 23 do begin  { 23 odpowiada –64 dB }
        db:=round(20*log10(i/35614));
        if db <> p then begin
          t[-1-db]:=i+1;
          p:=db
        end
      end;
  • #43
    slawek55
    Level 23  
    Dzięki.

    T aTwoja tabelka i funkcja decybel jest odwrotnie proporcjonalna do wartości wejsciowej liczbie S na wejsciu dużej odpowiada mała na wyjsciu funkcji.
    Jak robić tak aby liczbie np z końca zakresu (35614) odpowiadała liczba na końcu skali (64), bo narazie jest odwrotnie?

    A i skąd pomysł n takie wyliczenia z takelką - świetny!

    P.S.
    A jak powstała liczba 23 odpowiadająca -64db. U mnie ta liczba 64 to poprostu ilośc (rozdzielczośc) słupków.
  • #44
    Xitami
    Level 29  
    Decybel to stosunek. Możemy odnosić się do maksimum, wtedy 0 mam dla wartości maksymalnej, a minus ileś tam dla minimalnej.

    W funkcji decybel popraw na "return(64-j)" i będzie klawo. 64 dla max i 0 dla maleńkich.

    Napisałem, że sposób dobry do -77 dB, ale trzeba pamiętać jeszcze o czymś.
    Jeżeli przetwornik ma powiedzmy 10 bitów wtedy dynamika to tylko około 60 dB (bitowość*6).
    Wszystko co mniejsze to szum kwantyzacji, błędy zaokrągleń i promieniowanie reliktowe.

    20*log(23/35614) == -64 dB
    No to mamy 64 ziarenka w słupku co 1 dB.

    A pomysł z tabelką? Nie jest mój, prawie zawsze gdy chcemy szybko liczyć cos złożonego,
    można użyć tabelki z przygotowanymi wcześniej obliczeniami, prawie bo tabelka może być
    zbyt wielka, nie dla RAM ale dla CACHE (co nie grozi AVRom) i czasem na PC np. w FFT warto
    liczyć sinusa i cosinusa
  • #45
    r06ert
    Level 25  
    Witam! Sam postanowiłem się pobawić DFT ...i w końcu trafiłem tutaj ;)
    Napisałem taki fragment kodu:
    Code:

    #include <avr/io.h>
    #include <avr/pgmspace.h>
    #include "ffft.h"

    void fft_input (const int16_t *, complex_t *);
    void fft_execute (complex_t *);
    void fft_output (complex_t *, uint16_t *);
    //int16_t fmuls_f (int16_t, int16_t);

    int16_t capture[FFT_N];
    complex_t bfly_buff[FFT_N];   
    uint16_t spectrum[FFT_N/2];   

    int maina(){
    while(1){
    //pobieranie próbek

    fft_input(capture,bfly_buff);
    fft_execute(bfly_buff);
    fft_outpu(bfly_buff,spectrum);

    //wizualizacja wyniku
    }
    }


    jednak wyskakuje błędąd kompilacji:
    ffft.h:15: - error: multiple storage classes in declaration of `complex_t'

    już nie mam pomysły co to może byc... Czy któryś z szanownych kolegów którym udało się użyć biblioteki ffft.h pomógłby mi? z góry dzięki, pozdrawiam!
  • #46
    long_os
    Level 9  
    A mi podczas kompilacji Avr-gcc wypluł coś takiego:

    Code:
    > "make.exe" all
    

    -------- begin --------
    avr-gcc (WinAVR 20100110) 4.3.3
    Copyright (C) 2008 Free Software Foundation, Inc.
    This is free software; see the source for copying conditions.  There is NO
    warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


    Compiling C: main.c
    avr-gcc -c -mmcu=atmega8 -I. -gdwarf-2 -DF_CPU=16000000UL -Os -funsigned-char -funsigned-bitfields -fpack-struct -fshort-enums -Wall -Wstrict-prototypes -Wa,-adhlns=./main.lst  -std=gnu99 -MMD -MP -MF .dep/main.o.d main.c -o main.o
    In file included from main.c:33:
    ffft.h:12: error: multiple storage classes in declaration specifiers
    ffft.h:21: error: expected declaration specifiers or '...' before 'complex_t'
    ffft.h:24: error: expected ')' before '*' token
    ffft.h:25: error: expected ')' before '*' token
    main.c:37: error: redefinition of 'struct _tag_complex_t'
    main.c:40: error: 'complex_t' redeclared as different kind of symbol
    ffft.h:15: error: previous declaration of 'complex_t' was here
    main.c:42: error: conflicting types for 'fft_input'
    ffft.h:21: error: previous declaration of 'fft_input' was here
    make.exe: *** [main.o] Error 1

    > Process Exit Code: 2
    > Time Taken: 00:02


    Kod bez żadnych zmian i przeróbek :/ Co jest przyczyną ?