Elektroda.pl
Elektroda.pl
X
Proszę, dodaj wyjątek www.elektroda.pl do Adblock.
Dzięki temu, że oglądasz reklamy, wspierasz portal i użytkowników.

Potrzebna pomoc w Dyskretnej Transformacie Fouriera

slawek55 11 Wrz 2006 06:19 15853 45
  • #1 11 Wrz 2006 06:19
    slawek55
    Poziom 23  

    cześć.
    Mam prośbę czy może mi ktoś pomóc zrozymieć zastosowanie algorytmu DFT w uC AVR np Atmega16. Potrzebna mi pomoc w jego zastosowaniu w uC. Mam parę pytań.

    0 29
  • #2 11 Wrz 2006 10:46
    Zaquadnik
    Poziom 27  

    Sprecyzuj te pytania ;) Poza tym, o wiele lepsza i wydajniejsza obliczeniowo jest FFT. Taka odmiana DFT :]

    0
  • #3 11 Wrz 2006 10:59
    slawek55
    Poziom 23  

    Chcę zrobić analizator widma na procesorze Atmega32.
    użyję przetwornika A/c wbudowanego w niego i będę miał tablicę próbek sygnałów o np. N=512 elementach. Chcę poddać teraz ten zestaw próbek analizie DFT. MAm wzory matematyczne w których potrzebuję wyjaśnień.
    mogę liczyć na pomoc

    0
  • #4 11 Wrz 2006 21:27
    And!
    Admin grupy Projektowanie

    DFT może być prostsze w zaimplementowaniu na AVR niż FFT.
    DFT będzie łatwiej zrozumieć i jest dobrym wstępem do FFT.
    W EP leci teraz cykl o DFT (jak masz niedaleko bibliotekę to możesz się z tym zapoznać)

    Dość obszerny temat o FFT:
    https://www.elektroda.pl/rtvforum/viewtopic.php?t=14925
    oraz tutaj:
    https://www.elektroda.pl/rtvforum/topic182376.html

    DFT od deski do deski:
    http://vitecvitec.webpark.pl/dft/dft.htm
    strona główna:
    http://vitecvitec.webpark.pl/

    0
  • #5 12 Wrz 2006 09:11
    Xitami
    Poziom 29  

    Pewnie Sławku
    Choć 512 to dla maleńkiego avr'a trochę sporo
    256 w kilka milisekund? bardzo proszę w C i ASM

    0
  • #6 12 Wrz 2006 10:35
    slawek55
    Poziom 23  

    CZeść.
    Czy miało być coś załączone.
    CZy mogę się Ciebie troszke podpytać jak to zrealizować?

    0
  • #7 13 Wrz 2006 06:28
    slawek55
    Poziom 23  

    CZeść.
    Analizując opisy DFT mam wątpliwości:
    1. Jeśli pobieram sygnał z pzetwornika A/c to dlaczego trzeba każdą taką próbkę podzielić przez połowę rozdzielczości przetwornika. W opisie jest to określone jako odjęcie od wyniku składowej stałem. A dlaczego? CZy nie można liczyć na tych próbkach co są z A/C. Dlaczego się odejmuje wartośc stałej?
    2. Jaka jest wielkość bufora Re oraz Im i jak obliczyć amplitudy prążków widma. Raz znalazłem że Re podzielić przez N-ilość próbek wejściowych a raz że przez N/2

    A tak ogólnie to o co chodzi z tą składową stałą?

    0
  • #8 13 Wrz 2006 11:40
    slawek55
    Poziom 23  

    Cześć.
    Dalej próbuję walczyć z DFT.
    Dlaczego jak podaję się sygnał na przetwornik A/c w uC Atmega32 to od próbki odejmuje sie składową stałą. W książce Atmega w praktyce po odebraniu wynkiu z A/C odejmuje się od niegowartość 1024/2.
    Czy licząc DFT zawsze trzeba odejmować składową stałą. A tak w ogóle dlaczego się ją odejmuje.

    0
  • #9 13 Wrz 2006 17:57
    And!
    Admin grupy Projektowanie

    Po przeprowadzeniu analizy DFT w tablicy wyników pod indeksem 0 masz wartość składowej stałej. Na kolejnych pozycjach tablicy masz moce na kolejnych fragmentach widma próbkowanego sygnału.
    Składową stałą możesz się nie przejmować istotne są dane o zawartości częstotliwości składowych analizowanego sygnału.

    0
  • #10 13 Wrz 2006 18:39
    shg
    Specjalista techniki cyfrowej

    A masz gotowca :P

    Zaimplementowana jest FFT, długość bufora wynosi 128 próbek.
    Bufory są dwa, podobnie jak wejścia.
    Do ADC0 Doprowadza się sygnał przepuszczony przez filtr antyaliasingowy (dolnoprzepustowy) 16kHz (próbkowanie 32kHz, przetaktowany ADC), takie też jest pasmo analizatora.
    Do wejścia ADC1 doprowadza się sygnał przepuszczony przez filtr 2kHz (próbkowanie 4kHz).
    Bufor wejściowy wypełniany jest naprzemiennie, raz z jednego, raz z drugiego wejścia.
    Odświeżanie jest rzędu 20 razy na sekundę. Wąskim gardłem jest próbkowanie z kanału ADC1, gdyby dwukrotnie zwiększyć częstotliwość próbkowania na tym kanale to szybkość odświeżanie wzrośnie do około 30 razy/s.
    Większość analizatorów widma wyświetla prążki z FFT jak leci. Tu jest inaczej. Prążki są poddawane sumowaniu wagowemu, tak aby symulowane było zachowanie się analizatora analogowego (każdy z prążków jest odpowiednikiem sygnału przepuszczonego przez filtr Butterwortha 16 rzędu) Częstotliwości środkowe filtrów są niestandardowe, ale to tylko dlatego, że dostałem gotową matrycę LED do której musiałem sie dopasować. Bez problemu oczywiście możńa zrobić to zgodnie z normą. W plikach bars.h są wagi dla poszczególnych prążków, jest to po prostu macież (a właściwie to dwie, bo są dwa kanały) transmitancji filtrów analogowych dla każdego prążka otrzymanego z FFT. Wiersze odpowiadają prążkom wyświetlanym, kolumny to prążki z FFT.
    Zarówno skala częstotliwości, jak i amplitudy są logarytmiczne.
    Na skali Y można wybrać 6, lub 3 dB/punkt. Ponadto w ramach opcji jest też powolne opadanie prążków i markery wartości szczytowej (jak w analizatorze w winampie)
    Projekt jest jeszcze niedokończony, ale już niedługo. W zasadzie brakuje tylko klawiatury do sterowania wszystkim i tor analogowy trzeba nieco uprościć, bo obecnie filtry wejściowe wymagają symetrycznego zasilania. Prawdopodobnie całkowicie usunę filtr na wejściu ADC0, zostawie tylko jakiś prosty RC. Natomiast na wejściu ADC1 filtr o stromej charakterystyce jest koniecznością, tu zastosowałem filtr eliptyczny ósmego rzędu (MAX297), ale na wzmacniaczach operacyjnych też coś może być.
    Matryca podłączona jest przez dekoder BCD -> 1 z 10 (kolumny - dwie grupy po 6 kolumn przełączane tranzystorami), wiersze bezpośrednio do proca
    A tutaj można zobaczyć analizator w akcji, tyle że drochę się dźwięk z obrazem rozjechał przy montażu: http://shg.ovh.org/temp/fft/analizator.html

    Schemat. Niekompletny, bo jeszcze muszę tor analogowy przetestować: https://obrazki.elektroda.pl/27_1158165998.png

    Procedurę liczącą FFT znalazłem tu: http://elm-chan.org/works/akilcd/report_e.html

    0
  • #11 13 Wrz 2006 20:09
    slawek55
    Poziom 23  

    Rany jak ty do tego doszedłeś????

    Czy jest możliwość taka abyś opisał algortym FFT i dlaczego FFt a nie DFT

    0
  • #12 13 Wrz 2006 20:42
    Tdv
    Poziom 33  

    And! napisał:
    Po przeprowadzeniu analizy DFT w tablicy wyników pod indeksem 0 masz wartość składowej stałej. Na kolejnych pozycjach tablicy masz moce na kolejnych fragmentach widma próbkowanego sygnału.
    Składową stałą możesz się nie przejmować istotne są dane o zawartości częstotliwości składowych analizowanego sygnału.

    No z tym bym się nie zgodził. Odejmowanie 1024/2 to nie składowa stała sygnały, która faktycznie będzie w elemencie zerowym ale przesunięcie sygnału do pomiaru. Zwróć uwagę, że przetwornik w AVRach jest unipolarny, a przebiegi zmienne trzeba przetwarzać bipolarnym albo przesunąć sobie cały przebieg składową stałą do góry.
    Przykład (dla uproszczenia zakładam napięcie odniesienia 5V): sinusoida o amplitudzie 2V, normalnie próbka z 0 stopni ma 0, z 90 stopni ma 2, z 180 stopni ma 0 i z 270 stopni ma -2. -2 na przetworniku zabudowanym w AVR nijak nie osiągniesz, rozwiązanie dodajesz do sinusa napięcie stałe 2,5V i teraz probka dla zera ma wartość 2,5, dla 90 4,5, dla 180 ma 2, a dla 270 ma 0,5.
    Obróbka takich danych dała by błędny wynik, bo w zerowym elemencie byłoby 2,5, a mierzony sygnał ma składową stałą 0, te 2,5V to tylko taki trick do pomiaru i właśnie te 2,5V trzeba "cyfrowo" odjąć. W praktyce te 2,5V to dokładnei rozdzielczość przetwornika dzielona przez dwa.

    Podobny wątek już był:
    https://www.elektroda.pl/rtvforum/viewtopic.ph...32+widma&sid=1a15d36276c0543adb70a8cc413d4fb1

    0
  • #13 13 Wrz 2006 21:28
    shg
    Specjalista techniki cyfrowej

    Bo FFT jest duuuuużo szybsze.
    Złożoność obliczeniowa DFT to O(N^2), FFT to tylko N*log(N), gdzie N - ilość próbek. Wraz ze wzrostem N złożoność obliczeniowa FFT rośnie dużo wolniej niż bezpośredniej implementacji DFT.

    Dla przykładu bezpośrednia implementacja DFT dla 128 próbkek to 128^2 = 16384 operacje. FFT to około 270 operacji. Dla 1024 próbek to ospowiednie ponad milion i nieco ponad 3 tysiące. Jak widać im więcej próbek tym korzyści większe.
    Nie jest to do końca prawda, ale chodzi o samo porownanie i sens złożoności obliczeniowej.
    Ponadto mniejsza ilość operacji oznacza większą dokłądność w obliczeniach numerycznych (z powodu kumulowania się błędów wynikających ze skończonej precyzji obliczeń)

    FFT raczej nie opiszę, sam algorytm jest conajmniej bardzo skomplikowany, w ogole jest to jeden z najbardziej skomplikowanych algorytmów będących w powszechnym użyciu.

    Sczegółowy opis FFT można znaleźć np. w Eleanor Chu; Alan George, "Inside the FFT Blackbox, Serial and parallel Fast Fourier Transform algorithms"
    Od czytania tego może conajwyżej głowa rozboleć :P
    W zasadzie do wykorzystania FFT nie jest potrzebna znajomość zasady działania. Albo inaczej, zrozumieć jak działa FFT jest dużo łatwiej, niż zrozumieć dlaczego działa.
    Myślę, że można by zacząć stąd: http://en.wikipedia.org/wiki/FFT

    W telegraficznym skrócie FFT działa tak, że dzieli wejściowy bufor na mniejsze fragmenty.
    Konkretniej w algorytmie DIT (Decimation In Time) wygląda to tak, że wejściowy bufor dzielony jest na dwa dwuktornie mniejsze.
    Do pierwszego trafiają próbki nieparzyste, do drugiego parzyste.
    Potem każdy z tych buforów jest ponownie w ten sam sposób dzielony na dwie części i taki podział wykonuje się aż zostaną same wektory zawierające po dwie próbki (wtedy mówi się że radix = 2) i na nich wykonuje DFT, która dla dwóch próbek jest banalna (w obliczeniach bierze się wtedy sinus zera i PI radianów, mnożenie przez zera i jedynki). Podobnie można postąpić dla 3, 4 i 5 próbek (radix 3, 4 i 5), bo wartości sinusów i cosinusów dla kątów będących wielokrotnością 1/3, 1/4 i 1/5 2PI radianów można sobie wyliczyć, a znowóż dla radix = 4 są to proste (z punktu widzenia obliczeń) wartości).
    Właściwie to FFT z założenia dzieli bufor aż zostaną wektory jednoelementowe, a DFT pojedynczej wartości jest właśnie tą wartością. Tyle że ze względu na wydajność zazwyczaj podział kończy się na wektorze o długości 2 lub więcej.
    Wykorzystując transformaty z radixem niebędącym potęgą dwójki można wykonywać FFT na buforach o długości różnej od potęgi 2. Można nawet łączyć transformaty o różnych radiksach. Jest taka dość przyjazna implementacja, nazywa sie kissfft (do znalezienia na sourceforge), która działa właśnie na takiej zasadzie, dobiera różne wartości radiksa rozkładając zadaną długość na kolejne czynniki tak, aby były one jak największe (wybierane sa zpośród dostępnych) na przykłąd dla N=128 kolejne radiksy będą wynosić: 4, 4, 4 i 2 (bo 4*4*4*2 = 128), dla N = 100 będą to: 5, 5, 4 (bo 5*5*4 = 100)




    Ale najbardziej popularny jest radix = 2, ze względu na prostotę implementacji i niewielką objętość programu.
    Potem widma (czyli wyniki DFT) takich kawałków 'składa' się razem mnożąc przez tzw. twiddle factors, czyli zespolone pierwiastki n-tego stopnia z jedynki . Z punktu widzenia obliczeń są to po prostu sinusy i cosinusy odpowiednich kątów (w artykule na wikiepdii są linki, które to opisuja).
    Oprócz DIT są też algorytmy DIF (Decimation In Frequency), które dokonują zamiany miejscami próbek w wynikowym widmie. Ta zamiana próbek, o której piszę to właśnie to rozdzielanie bufora na mniejsze części. Najczęściej wykonuje się to tak że próbki zamienia się po prostu miejscami, tak że pierwsza połowa bufora to próbki nieparzyste, a druga - parzyste (i znowu to samo na pierwszej i drugiej połówce) bez fizycznego tworzenia kolejnych buforów. Proces ten nosi nazwę bit-reversal, ponieważ próbki zamieniane są miejscami w taki sposób, że gdyby ich indeksy zapisać w postaci binarnej, to okaże się że po operacji próbki znajdą się w takich pozycjach, jakby w ich indeksach cyfry zostały zapisane w odwrotnej kolejności. Na przykłąd próbka nr. 18 (binarnie: 00010010) w buforze o długości 256 otrzyma indeks 72 (binarnie 01001000).
    Ogólnie Algorytmy DIT i DIF nie różnią się specjalnie szybkością. Zresztą można nawet napisać algorytm, który będzie jednocześnie DIT i DIF, po prostu część operacji zostanie wykoanna w domenie czasu, a część w domenie częśtotliwości, ale dla przejrzystości raczej się tego nie praktykuje.
    Ponadto algorytmy FFT można podzielić na in-place i (jak w moim przypadku) out of place. Te pierwsze są nieco wolniejsze, ale wymagają mniej pamięci ponieważ widmo zostaje zapisane w miejscu gdzie wcześniej znajdowały się próbki sygnału w domenie czasu. W drugim wypadku widmo zapisywane jest do osobnego bufora i zawartość bufora z próbkami nie jest niszczona, można to też wykorzystać (tak jak u mnie) do przyspieszenia działania. Podczas gdy analizowany jest sygnał z pierwszego wejścia bufor próbek nie jest już potrzebny i jest wypełniany próbkami z drugiego wejścia, dzięki czemu program nie musi czekać bezproduktywnie na zakończenie obliczeń.
    Jeżeli natomiast masz procesor z dużą ilością rejestrów to wydajniejsze są algorytmy, które mają większy radix. Nie, AVR nie ma dużo rejestrów, bo raz - są one tylko ośmiobitowe, a dwa - nie wszystkie rejestry są równocenne, tzn. nie na każdych można wykonywać np. operacje mnożenia itp.
    Dla AVR najlepiej stosować radix = 2.

    Zresztą tak by można pisać i pisać. Pogrzeb w tych linkach z artykułu z wikipedii, znajdziesz tam między innymi obrazki, które dość sporo tłumaczą (na przykład "motylki" :]).

    W tej (bardzo dobrej zresztą) książce (można pobrać za free) też znajdziesz opis FFT, w miarę przyjazny: http://www.dspguide.com/

    A jak doszedłem? Posiedziałem, posiedziałem aż w końcu wydumałem. :P

    0
  • #14 14 Wrz 2006 07:10
    slawek55
    Poziom 23  

    Czym i jaka wersją kompilowałeś, bo mi w WinAVR wyskakują same błędy!! :-(

    Jeszcze jedno sobie przypomniałem. Odpowiedzieliście mi na pytanie o odjęciu warości1024/2 od wyniku a/c a dlaczego tu odejmujecie 32768 Skąd ta liczba, i co będzie jak sie nie odejmie tej liczby?

    0
  • #15 14 Wrz 2006 20:51
    shg
    Specjalista techniki cyfrowej

    slawek55 napisał:
    Czym i jaka wersją kompilowałeś, bo mi w WinAVR wyskakują same błędy!! :-(

    Najnowszą wersją WinAVR. W załączniku jeszcze raz to samo, ale kompletny projekt z Programmers Notepad. Przed chwilą sprawdziłem i kompiluje się tak jak powinno ('Make All' z menu Tools). Te procedury FFT mogą stwarzać problemy przy kompilacji, jak próbowałem ich pierwszy raz użyć, to się niesamowite cuda działy (niby wszystko wyglądało dobrze, ale tak jakby kompilator zapominał o deklaracjach).

    slawek55 napisał:
    Jeszcze jedno sobie przypomniałem. Odpowiedzieliście mi na pytanie o odjęciu warości1024/2 od wyniku a/c a dlaczego tu odejmujecie 32768 Skąd ta liczba, i co będzie jak sie nie odejmie tej liczby?

    Bo wynik z przetwornika A/C jest przesunięty w lewo, do 'krawędzi' (bit ADLAR ustawiony), czyli pomnożony przez 64.

    Wkleję tu, co dostałem na PW, bo może jesze komuś się przyda przy uruchamianiu.
    slawek55 napisał:
    Pomyślałem że moge podejść do sprawy inaczej.
    Skoro na tej stronie www http://elm-chan.org/works/akilcd/report_e.html, są gotowe funkcje do FFT to po co pisać je od nowa skoro można z nich skorzystać. Nie wiem czy zrobiłeś tak, ale dla mnie jest to jakieś rozwiązanie.

    Dokładnie tak zrobiłem. Pisanie własnej FFT to w zasadzie prawdziwe wyzwanie, wiele podręczników zresztą zaleca właśnie takie podejście, jeżeli jest coś gotowego, co spełni Twoje oczekiwania, to w zasadzie nie ma co kombinować. Własną FFT pisze się raczej tylko jak potrzebne są jakieś specyficzne jej własności. Na początku tak włąśnie chciałem zrobić - napisać FFT in-place dla wartości rzeczywistych (FFT jest zdefiniowana dla liczb zespolonych, dla liczb rzeczywistych można algorytm 'okroić' uzyskując przy tym około dwukrotne przyspieszenie), a in-place ze względu na ilość zajmowanej pamięci. Potem nawet znalazłem gotowy algorytm (algorytm Sorensena), tyle że na liczbach zmiennoprzecinkowych, których wykorzystanie nie wchodziło w grę ze względu na wydajność. W celu przystosowania praktycznie każdego algorytmu FFT do obliczeń na liczbach stałoprzecinkowych (takie już bez problemu 'strawi' AVR) żeby nie było strat na precyzji obliczeń należy wszystkie obliczenia ('wewnątrz') przeprowadzać z odpowiednim skalowaniem. W implementacji zmiennoprzecinkowej takiego problemu nie a, ponieważ można sobie przeskalować dopiero wynik końcowy.
    Niestety w wypadku wyżej wymienionego algorytmu okazało się to skrajnie trudne, jeśli nie niewykonalne. Wszystko przez nieco odmienny układ 'motylków' niż w bardziej 'standardowych' algorytmach, gdzie po prostu skalowanie wykonuje się przy danej operacji za każdym razem. W algorytmie Sorensena trzeba by wprowadzić dość skomplikowane reguły skalowania.

    slawek55 napisał:
    Moj problem teraz to jak to właśnie wykorzystać.
    S tam zarówno funkcje fft jak i do obsługi transmisji szeregowej, ale to akutat nie jest potrzebne.
    Czy możesz mi pomóc jak mogę to wykorzystć, bo prowdę mówiąc nie wiem co tam mam do wykorzystania i jak z tego skorzystać. Jaki rezonator tam jest gdzie wejście (tablica próbek) gdzie wyjście. Jakaś konfiguracja A/c itp.

    A no właśnie o tym nie pisałem. Częstotliwość rezonatora zdefiniowana jest w pliku Makefile (F_CPU) program prawie automatycznie dostosowuje się do tej częstotliwości (poza przetwornikiem A/C). Oczywiście im większa częstotliwość tym lepiej. W moim układzie jest 16MHz.

    Z w/w strony wykorzystałem tylko funkcje służące do liczenia FFT (fft_input(), fft_execute() i fft_output()) tabele motylków i okienka, oraz inne z nimi związane.

    Wszystkie peryferia inicjalizowane są w funkcji init().
    Co 31.25us wywoływane jest przerwanie timera 1 (OCR1A). Ten kawałek kodu wylicza wartość, którą trzeba wpisać do rejestru OCR1A timera
    Code:
    #define F_SMP      (32000UL)      /* sampling frequency */
    
    #define INT_CYCLES   (F_CPU / F_SMP)   /* interrupt period in timer cycles */

    Jeżeli aktualnie próbkowanie jest z wejścia ADC0 (pasmo 16kHz) to przy każdym przerwaniu pobierana jest próbka. Jeżeli próbkowanie jest z wejścia ADC1, to próbka pobierana jest co ósme przerwanie.
    Jeżeli bufor zostanie wypełniony to zerowana jest flaga BUFF_FILL w zmiennej flags.

    Teraz konfiguracja przetwornika. Ustawiony jest na tryb pojedynczej konwersji. Za każdym razem, gdy pobierana jest próbka (przerwanie timera 1, a konkretniej to OCR1A). Preskaler przetwornika ustawiony jest na 32 (dzieli zegar systemowy 16MHz przez 32 i w efekcie przetwornik taktowany jest częstotliwością 500kHz. Przetwornik jest przetaktowany (ale może poprawnie pracować z zegarem aż do 1MHz chyba), co zwiększa poziom szumu konwersji, no ale zwiększa też szerokość pasma. To jest ten fragment, który nie działa automatycznie (aczkolwiek nic nie stoi na przeszkodzie, poza chyba tylko moim lenistwem, żeby było inaczej), wartość preskalera należy dostosować do częstotliwości zegara samemu, tak aby zegar dla przetwornika wynosił minimum częstotliwość próbkowania * 13 (w tym wypadku (próbkowanie 32kHz) minimum to 416kHz)
    Code:
    /*
    
    ** enable ADC, prescaler = 16 -> ADC clock = 16MHz / 32 = 500kHz
    ** ADC is overclocked
    */
       ADCSRA = _BV(ADEN) + 5;


    Na początek możesz uruchomić samą FFT, nie jest to zbyt skomplikowane i ogólnie wygląda tak:

    Code:
    typedef struct _tag_complex_t {
    
       int16_t r;
       int16_t i;
    } complex_t;

    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];   /* Wave capture buffer */
    complex_t bfly_buff[FFT_N];         /* FFT buffer */
    uint16_t spectrum[FFT_N/2];      /* Spectrum output buffer */


    Pierwszy bufor (capture) to ten, do którego pobierane są próbki z przetwornika w ilości FFT_N.
    Drugi (bfly_buff) to bufor na którym wykonywana jest właściwa transformacja, dane z bufora capture trafiają tutaj po przepuszczeniu ich przez okienko Hamminga, dzięki temu częstotliwości obecne w sygnale, które nie trafiają dokładnie w częstotliwości otrzymywane w wyniku transformaty nie powodują powstawania 'rozlazłych' prążków (spectrum-leaks). Po okienkowaniu dane przenoszone są w inne miejsce (do bufora bfly_buff), więc bufor capture zostaje zwolniony i można rozpocząć próbkowanie z drugiego wejścia.
    Trzeci bufor to ten, do którego trafia widmo amplitudowe sygnału. Informacja o fazie jest odrzucana, ponieważ w analizatorze widma nie jest ona do niczego potrzebna. Stąd jest on o połowę krótszy.
    Algorytm bit-reversal aplikowany jest (w procedurze fft_output()), po wykonianiu FFT, a właściwie to tylko motylków (w fft_execute()), więc ta FFT jest algorytmem typu DIF. Tak właściwie to same motylki w tej implementacji wykonywane są in-place. Dopiero algorytm bit-reversal powoduje przeniesienie widma do kolejnego bufora.

    Dalej. Przetwornik A/C skonfigurowany tak, że wynik konwersji przesunięty jest w lewo (najstarszy bit wyniku konwersji to bit 15 w rejestrze ADC, lub 7 w ADCH, jak kto woli). Od tej wartości odejmowana (lub dodawana do niej, bez różnicy) jest wartość 0x8000, co powoduje konwersje z typu bez znaku na typ ze znakiem.
    W każdym przerwaniu po odczytaniu próbki z rejestru ADC natychmiast uruchamiana jest kolejna konwersja:
    Code:
       adc_smp = ((int16_t) ADC) - 0x8000;   /* read result of previous conversion */
    
       ADCSRA |= _BV(ADSC);               /* start new conversion */


    Gdy bufor zostanie zapełniony jego zawartość traktowana jest funkcją fft_input(), która aplikuje okienko i przenosi do kolejnego bufora - bfly_buff.

    Na tym buforze wykonywane są motylki (funkcja fft_execute())

    A potem bit-reversal i liczona wartość amplitudy/mocy, czy cokolwiek to jest (magnitude), jest to po prostu moduł wartości każdego prążka (DFT, a więci FFT daje widmo w postaci liczb zespolonych)


    W najprostszej wersji wygląda to tak:
    Code:
          fft_input(capture, bfly_buff);
    
          fft_execute(bfly_buff);
          fft_output(bfly_buff, spectrum);


    W spectrum dostaniesz widmo z liniową skalą częstotliwości i amplitudy.

    Funkcja pack_spectrum() dokonuje konwersji rozkładu częstotliwości z liniowego na logarytmiczny przez dosłownie 'upakowanie' pewnej ilości prążków otrzymanych z FFT do jednego prążka, który będzie wyświetlany. Na obrazku masz charakterystyki filtrów analogowych, które symyluje ten analizator.
    https://obrazki.elektroda.pl/62_1158259135.png
    Oś X to częstotliwość znormalizowana do częstotliwości próbkowania (czyli 1 na wykresie odpowiada częstotliwości 32kHz). Oś Y to wzmocnienie filtru dla danej częstotliwości.
    W programie wartości z wykresów zapisane są w tablicy, dla każdej częstotliwości prążka z FFT odczytywana jest wartość wzmocnienia danego filtru, amplituda prążka mnożona jest przez to wzmocnienie. Dla każdego filtru wartości te są sumowane tworząc w ten sposób wartość amplitudy, jaką uzyskalibyśmy na wyjściu filtru analogowego.
    Proces ten powtarza się dla każdego z filtrów.
    W bardziej matematycznym ujęciu wygląda to tak, że wektor amplitudy otrzymany z FFT mnożony jest przez macież transmitancji filtrów analogowych dając w efekcie wektor amplitudy odpowiadający filtrom analogowym.

    Z tak przekształconego widma funkcja make_image() tworzy obraz, który będzie wyświetlony na matrycy LED.
    W pierwszej kolejności Skala amplitudy widma konwertowana jest z postaci liniowej na logarytmiczną.
    Code:
          bar_height = 16 - clz(sig) + WINDOW_GAIN;

    sig to amplituda w postaci liniowej. bar_height to wysokość paska na matrycy, już w postaci logarytmicznej (od 0 do 16).
    Logarytm liczony jest sprytnie, za pomocą funkcji clz() Funkcja ta liczy ilość zer nieznaczących (z przodu) w zapisie binarnym wartości. Sprytność polega na tym, że można to zrealizować bardzo szybko, a po odjęciu wyniku od ilości bitów przez jaką była reprezentowana wartość w postaci liniowej otrzymuje się logarytm o podstawie 2, co przy okazji znakomicie nadaje się do zastosowań audio, ponieważ pozwala na uzyskanie bardzo przyjaznego kroku 3, lub 6 dB na punkt.
    WINDOW_GAIN to 'wzmocnienie' wnoszone przez okienkowanie, tu jest równe jeden, jednostką jest -6dB (minus), dlatego wartość jest dodatnia, okienko Hamminga tłumi o około 6dB, podobnie jak wiele innych (poza flat-top). Tłumienie należy potem dodać do otrzymanej wysokości prążka żeby uzyskać pełną skalę (bez tego nie zapalały by się najwyższe punkty na matrycy).
    Dalej jest już tworzony tylko obraz, który będzie wyświetlany na matrycy (tablica leds).
    To co tam wstawisz zależy już od konkretnego hardware'u. Podstawą jest wartość bar_height, którą należy w jakiś sposób wyświetlić.

    slawek55 napisał:
    A tak przy okazji dlaczego masz u siebie zastosowane dwa przetworniki i do czego to służy.

    No więc tak. Wykonanie transformacji Fouriera na 128 próbkach daje nam 64 prążki widma. Prążki rozłożone są liniowo w domenie częstotliwości (pomiędzy dowolnymi dwoma sąsiednimi prążkami jest taka sama różnica częstotliwości). W związku z tym dla częstotliwości próbkowania 32kHz odstępy między prążkami będą wynosić 500Hz. Prążki są rozłożone od 0 do 16kHz, czo wynika z twierdzenia Shannona o próbkowaniu (w trochę zmienionej formie: Widmo sygnału próbkowanego z częstotliwością x Hz ma szerokość x/2 Hz).
    Ludzki słuch (podobnie jak wszystkie inne zmysły) ma charakterystykę logarytmiczną, a precyzyjniej to tak: poziom odczuwanego bodźca jest proporcjonalny do logarytmu z jego intensywności. Dotyczy to także częstotliwości, najprościej mówiąc, im mniejsza częstotliwość, tym rozdzielczość większa (tak somo odczyjemy zmianę 50->51Hz, jak 5000->5100Hz, mimo że w pierwszym wypadku różnica wynosi 1, a w drugim 100Hz). Podobnie w muzyce, najwięcej 'dzieje się' w zakresie niższych częstotliwości (100Hz - 2 kHz) Dlatego ten zakres zajmuje ponad połowę matrycy.
    Zaistniała więc konieczność zwiększenia rozdzielczości w zakresie niźszych częstotliwości. Jednym ze sposobów jest właśnie pobranie próbek z mniejszą częstotliwością próbkowania, dzięki czemu widmo ulega zawężeniu, ale i to samo dotyczy rozdzielczości - odstępy między prążkami zmniejszają się tyle razy, ile razy zmniejszy się częstotliwość próbkowania.
    Dlatego podczas próbkowania z drugiego wejścia (ADC1) częstotliwość próbkowania wynosi 4kHz, a rozdzielczość widma 2kHz/64 = 31.25Hz.
    Na podstawie danych z wejścia ADC1 są 'składane' prążki od 62 do 1.2kHz. Na podstawie danych z ADC0 prążki od 2.1kHz do 16kHz.

    Częstotliwości prążków wyświetlanych na matrycy to kolejno:
    62Hz, 103Hz, 171Hz, 284Hz, 470Hz, 777Hz, 1.2kHz, 2.1kHz, 3.5kHz, 5.8kHz, 9.8kHz i 16kHz.
    Jest to rozkład logarytmiczny, stosunek dwóch sąsiednich częstliwości jest zawsze stały (103/62 = 171/103 = 284/171 itd.), w przeciwieństwie do wyniku otrzymywanego z FFT, gdzie różnica między sąsiednimi prążkami jest stała.
    Według normy ISO-IEC ileśtam częstotliwości ustala się tak, że bierze się częstotliwość bazową 1kHz i kolejne częstotliwości są większe i mniejsze od poprzedniej o 1/1, 1/2, 1/3 itd. (zależy od tego ile prążków chcemy) oktawy. Tu nie miałem takiej możliwości, trzeba było wykorzystać wszystkie słupki na matrycy stąd takie częstotliwości dziwne. W zasadzie najważniejsze żeby tylko rozkład logarytmiczny została zachowany.

    Inny sposób na uzyskanie większej rozdzielczości to zwiększenie ilości próbek, a co za tym idzie rozmiaru transformaty, jednak w AVR z racji małej ilości pamięci i małej mocy obliczeniowej jest to niepraktyczne.

    Jeszcze taka dość ważna uwaga. W funkcji clz() znajduje się taki fragment kodu:
    Code:
       x = (uint16_t) (((uint32_t) x * x) >> 16);

    Podnosi on przekazaną jako argument wartość do kwadratu i dzieli przez 2^(ilość bitów)
    w ten sposób uzyskuje się rozdzielczość amplitudy widma 3dB/punkt.
    Jeżeli fragment ten zostanie usunięty, to rozdzielczość widma wynosi 6dB/punkt.

    Warto poeksperymentować, obie wartości dają interesujące wyniki. Przy 3dB/punkt słupki na matrycy podskakują nieco 'żywiej' niż przy 6dB/punkt :] Docolowo będzie możliwość wyboru za pomocą klawiatury.

    Funkcja ilog2() jest narazie zbędna i docelowo zastąpi funkcję clz(). Efekt w zasadzie prawie ten sam, ale będzie zwracała odrazu wartość logarytmu.

    ----------
    Update.
    Dorzuciłem jeszcze kilka filmów. Na tych nowych jest skala 3dB/punkt, włączone markery wartości szczytowych i powolne opadanie słupków.
    Na filmie prazki.mpg po kolei wszystkie częstotliwości środkowe filtrów analogowych. Przy wyższych częstotliwościach pojawiają się krzaczki, zwalę winę na karte dźwiękową (badziewna). Wyższych częstotliwości nie słychać, to z kolei wina kamery.

    http://shg.ovh.org/temp/fft/analizator.html

    1
  • #16 15 Wrz 2006 07:17
    slawek55
    Poziom 23  

    Super, dziękuję.

    Ten kod przeszedł bez problemu bez ostrzeże i błędów, jest super.

    Z tego co napisałeś to to odjęcie 0x8000 powoduje przesunięcie wyniku zm przetwornika, a co ze składową stałą, poprzednio to co ja się pytałem o odjęciu 1024/2 to teraz jest jakby załatwione poprzez liczbę 0x8000. Czyli ta operacja powoduje zarówno przesunięcie poziomu oraz uporządkowanie odczytanej wartości.
    A jak bym zastosował Lcd geraficzny 128x64 ze sterownikiem KS0108 jak na podawanej wcześniej stronie www.apetech.de.. Problem w tym że na tej stronie nie jest podane z jaką czestotliwościa procuje ATmega128 obsługujący lcd, brak schematu.

    A wracając do opisu
    jaki zakres mają wartości w tabeli spectrum, czy są to już wartości RMS czy amplitudy, tzn, czy trzeba je jezcze oblicać czy można wyświetlić.
    Poniważ piszesz

    Cytat:
    A potem bit-reversal i liczona wartość amplitudy/mocy, czy cokolwiek to jest (magnitude), jest to po prostu moduł wartości każdego prążka (DFT, a więci FFT daje widmo w postaci liczb zespolonych

    Więc robi to któraś z funkcji z pliku FFT.S czy to Ty musiśz to wykonać.

    A co robi funkcja fmuls_f.

    A czy dobrze rozumiem że w tablicy spetrum są wartości dla konkretnej częstotliwości i aby była reakcja taka jaka jest jak piszesz w filtrach analogowych to trzeba kilka połączyć w jeden.

    Na koniec wspomniałem o wyświetlaczu graficznym poniweważ patrząc na film z tej strony http://www.apetech.de/article.php?artId=7&nnId=17 daje to fajny efekt, może nie widoczny zbyt wieczorem ale bardzo przypomina analizator w Winamie.

    Pa i jeszcze raz prosze o wyrozumiałośc w mioch wypowiedziach. Układ projekt jest super ale muszę przyzanć że zagadnienie bardzo ciążkie, a zależy mi na osiągnięciu celu.

    Sławek.


    POmyślałem o takiej jescze rzeczy,, jeśli by zmienić na początku organizacją bitów aby wynik z przetwornika był przesuniąty do prawej strony to co wtedy.

    0
  • #17 15 Wrz 2006 16:00
    shg
    Specjalista techniki cyfrowej

    slawek55 napisał:
    Z tego co napisałeś to to odjęcie 0x8000 powoduje przesunięcie wyniku zm przetwornika, a co ze składową stałą, poprzednio to co ja się pytałem o odjęciu 1024/2 to teraz jest jakby załatwione poprzez liczbę 0x8000. Czyli ta operacja powoduje zarówno przesunięcie poziomu oraz uporządkowanie odczytanej wartości.

    Dokładnie tak Skadowa stała wynosi 1/2 napięcia odniesienia, nie ma tego na schemacie. jeden rezystor 100k z wejścia przetwornika do AREF, drugi do masy, wejście sygnału przez kondensator.

    slawek55 napisał:
    A jak bym zastosował Lcd geraficzny 128x64 ze sterownikiem KS0108 jak na podawanej wcześniej stronie www.apetech.de.. Problem w tym że na tej stronie nie jest podane z jaką czestotliwościa procuje ATmega128 obsługujący lcd, brak schematu.

    http://www.apetech.de./article.php?artId=1&nnId=10
    Chyba o to chodzi?
    Jest tam taka tabelka:
    "Funtkion | Takte | Zeit bei 8MHz"
    Autor podaje czasy wykonywania funkcji przy zegarze 8MHz, więc chyba taki zegar należy zastosować (oczywiście im większy tym lepiej, ale trzeba by pewnie modyfikować nieco soft).

    slawek55 napisał:
    jaki zakres mają wartości w tabeli spectrum, czy są to już wartości RMS czy amplitudy, tzn, czy trzeba je jezcze oblicać czy można wyświetlić.
    Poniważ piszesz
    Cytat:
    A potem bit-reversal i liczona wartość amplitudy/mocy, czy cokolwiek to jest (magnitude), jest to po prostu moduł wartości każdego prążka (DFT, a więci FFT daje widmo w postaci liczb zespolonych

    Więc robi to któraś z funkcji z pliku FFT.S czy to Ty musiśz to wykonać.

    To co dostajesz po wykonaniu funkcji fft_output() to właśnie moduły wartości (to co pisałem jest przez tą funkcję wykonywane).
    W zasadzie jest to amplituda.
    Ponieważ wszystkie prążki są sinusami, to przeliczenie tego na wartości RMS sprowadza się jedynie do podzielenia amplitudy każdego prążka przez pierwiastek z 2.

    slawek55 napisał:
    A co robi funkcja fmuls_f.

    Mnożenie ułamkowe. Tak samo jak instrukcja FMUL/FMULS, tylko że szesnastobitowe.

    slawek55 napisał:
    A czy dobrze rozumiem że w tablicy spetrum są wartości dla konkretnej częstotliwości i aby była reakcja taka jaka jest jak piszesz w filtrach analogowych to trzeba kilka połączyć w jeden.

    Tak, w spectrum[] dostajesz kolejne częstotliwości rozłożone liniowo od 0Hz (składowej stałej) do częstotliwości próbkowania / 2.


    slawek55 napisał:
    POmyślałem o takiej jescze rzeczy,, jeśli by zmienić na początku organizacją bitów aby wynik z przetwornika był przesuniąty do prawej strony to co wtedy.

    Zmniejszy się dokładność obliczeń i końcowy wynik FFT będzie też mniejszy.

    0
  • #18 16 Wrz 2006 15:13
    slawek55
    Poziom 23  

    Cześć.

    Xastanawiałem się czym róznią się w sensie używania w programie dwie deklaracje:

    Code:
    void fft_input(const comlex_t *, comlplex_t *);
    
    od
    Code:
    void fft_input(const int16_t *, comlplex_t *);


    dodoatkowo są tam takie definicje
    Code:
    #define INPUT_NOUSE

    Code:
    #define INPUT_IQ

    ale jedna z nich jest nie używana. CZm one się róznią i gdzie są w programie zdefiniowane. Te dwie powyższe definicje są w pliku fft.h, ale sa poprzedzone znakami // wiec musi byc któraś z nich definowana gdzieś indziej.

    Zapomniałem, A co było z pierwszą wersją że wyskakiwało tyle błędów kompilacji?

    0
  • #19 16 Wrz 2006 19:46
    shg
    Specjalista techniki cyfrowej

    Ta nieużywana to wersja transformaty operująca na liczbach zespolonych na wejściu. Sygnał audio ma z natury charakter rzeczywisty, więc można stosować 'okrojoną' (szybszą) wersję transformaty.

    A co było z pierwszą wersją, to nie wiem, obie są takiew same, do tej drugiej dorzuciłem tylko plik projektu z programmers notepad. Może coś nie tak z konfiguracją.

    0
  • #20 16 Wrz 2006 20:20
    slawek55
    Poziom 23  

    Ok. Rozumiem.

    Powiedz mi taką rzecz, jaki zakres wartości przyjmuje spectrum po obliczeniach od ilu do ilu. Są to zmienne unsigned, tak. A jak to się ma do zmiennej capture z przetwornika

    0
  • #21 17 Wrz 2006 16:23
    slawek55
    Poziom 23  

    Próbowałem w AvrStudio4 jak przebiegają obliczenia na wynikach otrzymanych z przetwornika A/C.
    Chodzi o ten wzór

    Code:
    adc_smp = ((int16_t) ADC) - 0x8000;


    Nijak sie mi to ma do tego co wydaje sie że powinno wyjśc przy normalnych obliczeniach robionych na piechotę czyli przesunięcie z powrotem bitów w prawo o 6 pozycji i odjęcie składowej stałej (połowy przetwornika czyli 512).

    Powiem że nic z tego nie rozumiem o co tu chodzi, czy AvrStudio cos namieszał czy ja niewidzę jakiegoś szczegółu który pozwoli to wyjaśnić.

    0
  • #22 20 Wrz 2006 18:19
    shg
    Specjalista techniki cyfrowej

    slawek55 napisał:
    jaki zakres wartości przyjmuje spectrum po obliczeniach od ilu do ilu. Są to zmienne unsigned, tak. A jak to się ma do zmiennej capture z przetwornika


    Echhh. A nie pamiętam dokładnie, jak zaczynałem pisać program, to dorzuciłem kawałek kodu, który pokazał tą wartość na LCD. Wydaje mi się, że było to 16383. Dolna granica to zero, bo jest to moduł wartości zespolonej (odpowiednik wartości bezwzględnej dla liczb rzeczywistych), więc jest zawsze dodatni.

    W tablicy capture dane przyjmują wartości od -32768 do +32767 czyli pełny zakres 16 bitów ze znakiem.

    slawek55 napisał:
    Nijak sie mi to ma do tego co wydaje sie że powinno wyjśc przy normalnych obliczeniach robionych na piechotę czyli przesunięcie z powrotem bitów w prawo o 6 pozycji i odjęcie składowej stałej (połowy przetwornika czyli 512).


    Nie za bardzo rozumiem :/
    Z tego kawałka kodu dostajesz wartość ze znakiem w zakresie -32768..+32767, dlatego że po konwersji wynik przesuwany jest w lewo. Najmniej znaczące bity są wypełniane zerami.

    Jeżeli w tak otrzymanym wyniku przesuniesz bity w prawo, to dostaniesz sieczkę.
    Należy wykonać przesunięcie arytmetyczne, czyli z zachowaniem znaku. Różnisię ono od przesunięcia logicznego tym, że od strony MSb wstawiane są nie zera, a kopiowany jest najstarszy bit, czyli dla wartości ujemnych najstarsze bity zostaną wypełnione jedynkami, a dla zera i dodatnich zerami. Kompilator C dobiera automatycznie rodzaj przesunięcia w zależności od typu użytych zmiennych.

    Jeżeli z kolei chcesz odejmować składową stałą bez przesuwania, to albo przy konfiguracji przetwornika A/C nie ustawiasz bitu za to odpowiedzialnego i odejmujesz 512:
    Code:
    adc_smp = ((int16_t) ADC) - 0x0200;

    albo (jeżeli tą samą wartość chcesz uzyskać z tego kawałka kodu, który wkleiłeś) wykonujesz TYLKO przesunięcie (arytmetyczne), bo składowa stała została już usunięta.

    0
  • #23 20 Wrz 2006 20:08
    slawek55
    Poziom 23  

    Cytat:
    slawek55 napisał:
    Nijak sie mi to ma do tego co wydaje sie że powinno wyjśc przy normalnych obliczeniach robionych na piechotę czyli przesunięcie z powrotem bitów w prawo o 6 pozycji i odjęcie składowej stałej (połowy przetwornika czyli 512).


    Nie za bardzo rozumiem :/

    To już się wyjaśniło. Na początku myślałem, że odjęcie liczby 0x8000 powoduje odjęcie składowej stałej i przesuwanie bitów, ale nie, więc wszystko jasne.


    przeprowadziłem symulacje w AvrStudio 4, mi wyszło że maksymalna wrtościa jaka może pojawić sie w zmiennej spectrum to 35314.

    W tablicy capture, to co do liczby -32768 się zgodzę ale do dodatnich to nie. mi wychodzi +32704, bo przecież bity najmniej znaczące nie są ustawione (przy max watości podanej na A/C 1 są tylko na 10 bitach)

    A czy masz jakiś sposób aby tą max wartośc w tablicy spectrum przedstawić na lcd o rozdzielczości od 0 do ok 60.

    Dodano po 18 [minuty]:

    A przepraszam być może jest to ważne ale ja symulowałem podanie składowej stałej na wejście i ta liczba dotyczyła indeksu 0 w tablicy capture. CZy oze miec to jakieś znaczenie w stosunku do reszty w tablicy, tzn czy inne (od 1 to końca) sa inczej wyliczane mogą być mniejsze?

    0
  • #24 20 Wrz 2006 22:00
    shg
    Specjalista techniki cyfrowej

    slawek55 napisał:
    przeprowadziłem symulacje w AvrStudio 4, mi wyszło że maksymalna wrtościa jaka może pojawić sie w zmiennej spectrum to 35314.

    A całkiem możliwe, jak pisałem - sprawdzałem to dość dawno.

    slawek55 napisał:
    W tablicy capture, to co do liczby -32768 się zgodzę ale do dodatnich to nie. mi wychodzi +32704, bo przecież bity najmniej znaczące nie są ustawione (przy max watości podanej na A/C 1 są tylko na 10 bitach)

    No tak. To był taki skrót myślowy. Gdyby podpiąć przetwornik 16 bitowy, to by było 32767 i z takimi wartościami może poradzić sobie oprogramowanie. Z wbudowanego ADC jest mniej; tak jak napisałeś.

    slawek55 napisał:
    A czy masz jakiś sposób aby tą max wartośc w tablicy spectrum przedstawić na lcd o rozdzielczości od 0 do ok 60.

    W skali logarytmicznej należało by podnieść wartość z tablicy spectrum do 4 potęgi, tak żeby zajmowała 64 bity i logarytmować, rozdzielczość będzie wynosić wtedy 1.5dB na punkt. Do logarytmowania conajmniej wskazane było by użycie nieco sprytniejszego algorytmu, ten który jest działa dobrze dla 8 i 16 bitów, im więcej bitów, tym wolniej.
    Ogólnie ten sprytniejszy algorytm (funkcja clz())polega na sprawdzeniu jak największej ilości początkowych zer za jednym zamachem, czyli zamiast sprawdzać bit po bicie zaczynamy sprawdzanie od całych bajtów z każdym zerowym bajtem dodając 8 do ilości zer i dopiero jak trafimy na niezerowy bajt zaczynamy sprawdzanie po jednym bicie.
    Można by też od biedy w skali liniowej wyświetlać, ale nie polecam, paskudnie to wygląda.

    slawek55 napisał:
    A przepraszam być może jest to ważne ale ja symulowałem podanie składowej stałej na wejście i ta liczba dotyczyła indeksu 0 w tablicy capture. CZy oze miec to jakieś znaczenie w stosunku do reszty w tablicy, tzn czy inne (od 1 to końca) sa inczej wyliczane mogą być mniejsze?

    Wszytkie elementy tej tablicy traktowane są jednakowo.
    Ale jeżeli chcesz uzyskać w wyniku symulacji jakieś konkretne widmo, to skłądową stałą trzeba wprowadzać do wszystkich elementów.
    W przeciwnym wypadku będziesz miał w pierwszej próbce deltę Diraca, która daje w wyniku płaskie widmo w domenie częstotliwości.

    I tak jeszcze przy okazji:
    W wyniku symulacji filtrów analogowych w takim układzie wysokości prążków reprezentują całkowitą energię w pasmie przenoszenia danego filtru (większa dla wyższych częstotliwości, bo i pasmo szersze). Stąd po przetworzeniu delty Diraca prążki odpowiadające wyższym częstotliwościom pokażą większą wartość. Jest to jak najbardziej celowe, ponieważ rozkład amplitudy w muzyce (i ogolnie w naturalnych dźwiękach) to 1/f (fale o wyższej częstotliwości mają wyższą energię przy tej samej amplitudzie). Odpowiada to charakterystyce szumu różowego.
    Gdy na wejście analizatora poda się szum różowy, to wszystkie prążki na matrycy będą miały taką samą wysokość.

    0
  • #25 21 Wrz 2006 07:10
    slawek55
    Poziom 23  

    A możesz mi powiedzieć jak to zrobić z tym logarytmowaniem, o którym mi napisałeś.

    0
  • #26 22 Wrz 2006 00:18
    shg
    Specjalista techniki cyfrowej

    W C to by tak wyglądało:

    Code:
    uint8_t clz(uint64_t x) {
    
    uint8_t n=0;
       if(x == 0) /* To niekoniecznie, troche trwa takie sprawdzanie, ale zeby sie tego pozbyc trzeba by zmodyfikowac petle ponizej, zeby sie nie 'wysypaly' */
          return 64;
       
       /* Najpierw po calym bajcie */
       while((x & 0xff00000000000000ULL) == 0) {   /* dopoki najstarszy bajt = 0 */
          n += 8;      /* dodac 8 zer do licznika */
          x <<= 8;   /* wartosc przesunac o 8 bitow w lewo */
       }
       /* A teraz po jednym bicie */
       while((x & 0x8000000000000000ULL) == 0) {   /* dopoki najstarszy bajt = 0 */
          n += 1;      /* dodac 1 zero do licznika */
          x <<= 1;   /* wartosc przesunac o 1 bit w lewo */
       }
          
       return n;
    }


    Nie radzę tego próbować. avr-gcc jest mało 'kumatym' kompilatorem i cos mi się zdaje, że nie poradzi sobie z optymalizacją.
    To samo napisane w asm będzie bardzo szybkie.
    Generalnie optymalizacja polega na dostępie do zmiennej za pomocą wskaźnika na bajt, można to i w C zrobić, pewnie nawet jako tako zadziała.

    W asm to bedzie najpierw taka pętelka, z użyciem instrukcji LD Rd, X+ (big-endian), lub LD Rd, X- (little-endian), potem TST Rd, jakies instrukcje warunkowe i dodawanie do licznika.
    W pierwszym niezerowym bajcie zera liczone za pomocą petli z LSL i sprawdzaniem flagi C.

    0
  • #27 24 Wrz 2006 09:17
    slawek55
    Poziom 23  

    Cześć.
    A jestes w stanie mi jakoś pomóc, z ym logarytmem, przyznam, zę ciężko i to robić w asemblerze.
    A przy okazji, próbuję wykombinować jak dodoać do tego opadanie słupków z tzw. Pick Holding. Jaki jest algorytm tego. Możesz coś podpowiedzieć.

    0
  • #28 24 Wrz 2006 18:28
    Xitami
    Poziom 29  

    S = RE*RE + IM*IM
    dB=3*log2(S) – 87

    Code:
     int log2(long x) {
    
       char r=0;
       if ( x & 0xffff0000 ) { x >>=16; r +=16; }
       if ( x & 0x0000ff00 ) { x >>= 8; r += 8; }
       if ( x & 0x000000f0 ) { x >>= 4; r += 4; }
       if ( x & 0x0000000c ) { x >>= 2; r += 2; }
       if ( x & 0x00000002 ) {          r += 1; }
       return r; }

    S ==_hh:_hl:_lh:_ll (cztery rejestrry)
    Code:
       clr   _bl   
    
       mov   _al,_hh
       or   _al,_hl
       breq   _msb1
          movw   _ll,_hl
          subi   _bl,low(-16)
       _msb1:
       tst   _lh
       breq   _msb2
          mov   _ll,_lh
          subi   _bl,low(-8)
       _msb2:
       mov   _al,_ll
       andi   _al,0xf0
       breq   _msb3
          swap   _ll
          subi   _bl,low(-4)
       _msb3:
       mov   _al,_ll
       andi   _al,0x0c
       breq   _msb4
          lsr   _ll
          lsr   _ll
          subi   _bl,low(-2)
       _msb4:
          andi   _ll,2
          breq   _msb5
          subi   _bl,low(-1)
       _msb5:
       ;//l=j*3-87;         
       mov   _l,_bl   ;
       add   _l,_l   ;
       add   _l,_bl   ;
       mov   _bh,_l
       subi   _l,low(87) 

    wynik w _l, co 3 dB
    -----------------------------------------------------------------------
    Warto sprawdzić jak zrealizowano biblioteczną funkcję lsqrt(), pierwiastek z 32 bitów. Widziałem to czasem bardzo źle zrobione (długo się liczy).
    Code:
           ;//fr[i]=lsqrt(ti:tr);
    
       clr   _qi
       clr   _qih
       movw   _qr,_qi
       movw   _wi,_qi
       movw   _wr,_qi
       movw   r30,_qi
       mov   _istep,_qi     
       ldi   r27,low(16)
       sqrtstart:   
          lsl   _tr      ;r:x <<= 2; wi:wr:ti:tr<<=2
          rol   _trh
          rol   _ti
          rol   _tih
          rol   _wr
          rol   _wrh
          rol   _wi
             
          lsl   _tr
          rol   _trh
          rol   _ti
          rol   _tih
          rol   _wr
          rol   _wrh
          rol   _wi

          lsl   r30      ; a<<=1
          rol   r31
          movw   _qi,r30    ;e=(a<<1)|1
          sec
          rol   _qi
          rol   _qih
          cp   _wr,_qi      ;r>=e?
          cpc   _wrh,_qih
          cpc   _wi,_istep
          brlo   sqrtend
             sub   _wr,r12   ;r-=e
             sbc   _wrh,r13
             sbc   _wi,r14   
             adiw   r30,1   ;a++
          sqrtend:
          dec   r27
          brne   sqrtstart 
       ldi   r27,high(_fr)
       mov   r26,_i
       lsl   r26
       st   x+,r31         ;fr[i]=swap(sqrt()); !!!!!!
       st   x,r30
    --------------------------------------------------
    Code:
       ;.def   _tr=   r2
    
       .def   _trh=   r3
       ;.def   _ti=   r4
       .def   _tih=   r5
       ;.def   _wr=   r6
       .def   _wrh=   r7
       ;.def   _wi=   r8
       .def   _wih=   r9
       ;.def   _qr=   r10
       .def   _qrh=   r11
       ;.def   _qi=   r12
       .def   _qih=   r13
       ;.def   _istep=   r14

       .def   _param=   r15

       .def   _i=   r16
       .def   _j=   r17
       .def   _l=   r18
       .def   _m=   r19
       .def   _bl=   r20
       .def   _bh=   r21

       .def   _al=   r22
       .def   _ah=   r23
       .def   _k=   r24
       .def   _zero=   r25
       .def   _ll=   r26
       .def   _lh=   r27
       ;r28
       ;r29
       .def   _hl=   r30
       .def   _hh=   r31

    0
  • #29 24 Wrz 2006 18:35
    slawek55
    Poziom 23  

    A mogę cię prosić o jakiś komentarz do tego. Bardzo prosze

    0
  • #30 24 Wrz 2006 20:54
    shg
    Specjalista techniki cyfrowej

    Co do peak - hold, ja to mam tak:
    Kazdy slupek ma swoj licznik opoznienia i zmienną przechowujaca wartosc szczytowa.

    Po obliczeniu wysokosci slupka program sprawdza, czy jest ona wyższa od wartości szczytowej.
    Jeżeli tak, to wysokość jest kopiowana do zmiennej przechowującej wartość szczytową i ustawiany jest licznik opoznienia na żądaną wartość, powiedzmy 500ms. Załóżmy teraz, że przerwanai uaktualniające wyświetlacz wywolywane sa co 10ms, wiec licznik powyzej bedzie ustawiany na 50.

    W każdym przerwaniu wykonuje się następujące czynności:
    - Jezeli dla słupka licznik opoznienia ma niezerowa wartosc, to jest zmniejszany o 1.
    - Jezeli licznik ma wartosc zerową (upłynął czas wyświetlania wartości szczytowej), to co określony czas zmniejszana jest wartość szczytowa, o ile jest większa od zera. (powolne opadanie markera wartości szczytowej)

    Podobnie można zrobić ze słupkami widma (powolne opadanie). Jeżeli obliczona wysokość słupka jest wyższa od aktualnej, to wysokość jest uaktualniana.
    W przerwaniu co określony interwał czasu wartość wysokości zmniejszana jest o 1 (o ile jest większa od zera).

    ------------------

    Jeszcze jeden kod logarytmowania, myślę, że gcc sobie z tym w miarę dobrze poradzi, a nawet jak nie, to sam kod jest bardzo podatny na optymalizację:

    Code:
    unsigned int intlog2(unsigned int x) {
    
    int log2=16;
     if((x & 0x00ff) == 0) {
       log2 -= 8;
       x <<= 8;
     }
     if((x & 0x0fff) == 0) {
       log2 -= 4;
       x <<= 4;
     }
     if((x & 0x3fff) == 0) {
       log2 -= 2;
       x <<= 2;
     }
     if((x & 0x7fff) == 0) {
       log2 -= 1;
       x <<= 1;
     }
     if(x == 0)
       return(0);

     return(log2);
    }

    0