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
#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)
/*
** 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:
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:
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:
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ą.
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:
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