FFT - cd....
Dziś będzie nieco bardziej teoretycznie i programistycznie.
Ostatnio pokazywałem na rysunku symbolicznie, że complex fft pozwala na wyliczenie zarówno dodatnich jak i ujemnych częstotliwości w widmie odbieranego sygnału.
Dla przypomnienia, jeśli dane z detektora Tayloe (wydzielającego sygnały I oraz Q) są próbkowane z częstotliwością f to w sygnałach I oraz Q pojawią się jedynie częstotliwości co najwyżej równe +-f/2. Walet działa z częstotliwością f = 24kHz. Complex fft potrafi rozróżnić te częstotliwości dlatego przy obrazowaniu sygnału zastosowałem skalę zarówno na plus jak i na minus.
[Obrazek: https://lx-net.pl/walet/image/fft_c01d.jpg]
Zeru odpowiada częstotliwość odbioru a skala rozciąga się od -10kHz do +10kHz choć teoretycznie można byłoby ją pociągnąć do +-12kHz. Ale ponieważ poziome rozmiary skali to nieco więcej niż 200 pikseli więc łatwo jest dopasować jedno do drugiego kiedy obrazuje się sygnały w zakresie +-10kHz.
Oś pionowa to w przybliżeniu skala S gdzie odległości między wartościami S są w przybliżeniu równe ale dodatkowo, im silniejszy sygnał tym "jaśniejszy" ma kolor.
Transformata fft wylicza tzw. bin-y oznaczające w skrócie oszacowania zawartości częstotliwości wynikające z kolejnego parametru tzw. fft_size. fft_size zwykle ma wartość równą potędze 2 i w przypadku zastosowanym w Walecie to liczba 256 i odpowiada ona za liczbę bin-ów wyliczonych w transformacie.
W transformacie typu real uzyskuje się w praktyce połowę z fft_size a więc 128 różnych bin-ów a w transformacie complex jest ich faktycznie tyle ile fft_size - oddzielnie 128 dla częstotliwości mniejszych i oddzielnie 128 dla częstotliwości większych od częstotliwości odbioru. Ponieważ w obu przypadkach w widmie sygnału pojawi się co najwyżej częstotliwość +-12kHz to na jeden bin przypadnie 12kHz/128=93.25Hz. W zakresie 10kHz na wykresie trzeba więc zastosować 10000/93.25 = ok. 107 binów po stronie dodatniej i po stronie ujemnej.
Z natury obliczeń transformaty wynika, że jeśli w widmie odbieranego sygnału będzie częstotliwość odpowiadająca dokładnie wielokrotności częstotliwości bin-a a więc 93.75Hz, 187.5Hz, 281.25Hz ...... itd. to wyliczenie będzie skupiać się na jednej wartości i powstanie jeden prążek na ekranie ale im bardziej będzie od niej oddalona tym więcej "energii" będzie rozproszone na sąsiednie bin-y a obraz na ekranie nie będzie pojedynczym prążkiem ale "choinką" prążków.
Oczywiście, im więcej bin-ów (czyli fft_size) tym lepiej dla dokładności pomiarów widma sygnału ale wiąże się to z większym nakładem obliczeniowym.
I tu pojawia się algorytm szybkiej, dyskretnej transformaty Fouriera opracowany przez panów Cooley i Tukey jeszcze w roku 1965 wymagający aby ilość danych (fft_size) miała właśnie rozmiar potęgi liczby 2. W skrócie pozwala to na zmniejszenie ilości obliczeń (mnożeń) z N^2 do wartości N*log2(N). W przypadku Waleta gdzie N=256 jest to zysk z N^2 = 65536 do 256*log2(256)=256*8=2048 czyli 32 razy mniej. Dla większych wartości N zysk jest jeszcze bardziej widoczny.
[Obrazek: https://lx-net.pl/walet/image/motylek1.jpg]
Algorytm ten zwany jest również motylkowym lub radix-2 ze względu na graficzne podobieństwo zobrazowania par obliczeń do skrzydeł motyla. Na każdym kroku algorytmu łączy się do obliczeń ze sobą pary wartości stąd wymaganie aby ilość danych wejściowych była potęgą liczby 2.
A teraz nareszcie w skrócie na temat samego algorytmu complex fft zastosowanego w kodzie Waleta.
[Obrazek: https://lx-net.pl/walet/image/fftc_dsp.jpg]
To co poprzednio oznakowałem jako blok cfft w rzeczywistości składa się z dwóch głównych części. Pierwsza, skojarzona z podzielonymi przez MultiOutput danymi I oraz Q zwana iqReceiver z klasy CallbackStream pozwala na odczyt próbki 2*256 danych z kanałów I oraz Q a drugą część to właściwa complex fft składająca się z dwóch procedur: dsps_fft2r_fc32 oraz dsps_bit_rev_fc32. Na końcu uzyskuje się macierze wartości rzeczywistych i urojonych dla każdego bin-a, które służą do wyliczenia jego amplitudy na zasadzie składania wektorów: globalFullAmplitudes[k] = sqrtf(real_part * real_part + imag_part * imag_part). MultiOutput zapewnia dostarczenie danych I oraz Q oraz uruchomienie funkcji obliczeń transformaty po skompletowaniu 2*256 danych poprzez procedurę writeDataCallback.
W pliku tekstowym wydzieliłem najważniejsze polecenia dotyczące complex fft
Analiza danych z tego pliku pewnie zainteresuje tylko programistów ale muszę przyznać, że uzyskanie tego kodu jest na granicy moich możliwości programistycznych, szczególnie przy obsłudze danych i procedur dotyczących dużych struktur oraz funkcji i metod Arduino Audio Tools. Posiłkowałem się przy tym AI, która oczywiście halucynowała proponując nieraz dziwne rozwiązania i przepraszając za błędy ale w końcu udało się uzyskać poprawny wynik.
Ostateczna procedura complex fft nie jest częścią Arduino Audio Tools ale stanowi natywne opracowanie producenta procesora stąd na początek ładowana jest biblioteka #include "esp_dsp.h" a w dalszej części uruchamiane jest sprzętowe wspomaganie fft.
W wyniku działania opisanych procedur powstaje tabela globalFullAmplitudes zawierająca 256 binów ułożonych według następującego porządku: bin 0 to tzw. składowa stała, która w tym przypadku jest pomijana, bin 1 to wartość dla 93.5Hz, bin 2 to wartość dla 187Hz ..... bin 127 to wartość dla 12000Hz, bin 128 to wartość dla -12kHz, ...... bin 254 to wartość dla -187.5Hz i bin 255 to wartość dla częstotliwości -93.75Hz. Jak widać dane w tabeli są dla częstotliwości dodatnich i ujemnych podane w odwróconym porządku co należy uwzględnić w graficznym obrazowaniu tych danych.
Udało mi się również ostatnio uzyskać obrazowanie sygnału z pomocą fft również w formie "wodospadu" ale o tym napiszę w kolejnym poście jak przygotuję zarówno rysunki, zdjęcia oraz nagrania i obiecuję, że nie będzie więcej tak mocno programistycznych postów, choć oczywiście postaram się coś wyjaśnić w razie pytań i potrzeby.
Szykuję też pliki Arduino i binarne dla Github ale również potrzebuję trochę czasu na uporządkowanie kodu przez jego prezentacją.
cdn. lj
Github
Dziś będzie nieco bardziej teoretycznie i programistycznie.
Ostatnio pokazywałem na rysunku symbolicznie, że complex fft pozwala na wyliczenie zarówno dodatnich jak i ujemnych częstotliwości w widmie odbieranego sygnału.
Dla przypomnienia, jeśli dane z detektora Tayloe (wydzielającego sygnały I oraz Q) są próbkowane z częstotliwością f to w sygnałach I oraz Q pojawią się jedynie częstotliwości co najwyżej równe +-f/2. Walet działa z częstotliwością f = 24kHz. Complex fft potrafi rozróżnić te częstotliwości dlatego przy obrazowaniu sygnału zastosowałem skalę zarówno na plus jak i na minus.
[Obrazek: https://lx-net.pl/walet/image/fft_c01d.jpg]
Zeru odpowiada częstotliwość odbioru a skala rozciąga się od -10kHz do +10kHz choć teoretycznie można byłoby ją pociągnąć do +-12kHz. Ale ponieważ poziome rozmiary skali to nieco więcej niż 200 pikseli więc łatwo jest dopasować jedno do drugiego kiedy obrazuje się sygnały w zakresie +-10kHz.
Oś pionowa to w przybliżeniu skala S gdzie odległości między wartościami S są w przybliżeniu równe ale dodatkowo, im silniejszy sygnał tym "jaśniejszy" ma kolor.
Transformata fft wylicza tzw. bin-y oznaczające w skrócie oszacowania zawartości częstotliwości wynikające z kolejnego parametru tzw. fft_size. fft_size zwykle ma wartość równą potędze 2 i w przypadku zastosowanym w Walecie to liczba 256 i odpowiada ona za liczbę bin-ów wyliczonych w transformacie.
W transformacie typu real uzyskuje się w praktyce połowę z fft_size a więc 128 różnych bin-ów a w transformacie complex jest ich faktycznie tyle ile fft_size - oddzielnie 128 dla częstotliwości mniejszych i oddzielnie 128 dla częstotliwości większych od częstotliwości odbioru. Ponieważ w obu przypadkach w widmie sygnału pojawi się co najwyżej częstotliwość +-12kHz to na jeden bin przypadnie 12kHz/128=93.25Hz. W zakresie 10kHz na wykresie trzeba więc zastosować 10000/93.25 = ok. 107 binów po stronie dodatniej i po stronie ujemnej.
Z natury obliczeń transformaty wynika, że jeśli w widmie odbieranego sygnału będzie częstotliwość odpowiadająca dokładnie wielokrotności częstotliwości bin-a a więc 93.75Hz, 187.5Hz, 281.25Hz ...... itd. to wyliczenie będzie skupiać się na jednej wartości i powstanie jeden prążek na ekranie ale im bardziej będzie od niej oddalona tym więcej "energii" będzie rozproszone na sąsiednie bin-y a obraz na ekranie nie będzie pojedynczym prążkiem ale "choinką" prążków.
Oczywiście, im więcej bin-ów (czyli fft_size) tym lepiej dla dokładności pomiarów widma sygnału ale wiąże się to z większym nakładem obliczeniowym.
I tu pojawia się algorytm szybkiej, dyskretnej transformaty Fouriera opracowany przez panów Cooley i Tukey jeszcze w roku 1965 wymagający aby ilość danych (fft_size) miała właśnie rozmiar potęgi liczby 2. W skrócie pozwala to na zmniejszenie ilości obliczeń (mnożeń) z N^2 do wartości N*log2(N). W przypadku Waleta gdzie N=256 jest to zysk z N^2 = 65536 do 256*log2(256)=256*8=2048 czyli 32 razy mniej. Dla większych wartości N zysk jest jeszcze bardziej widoczny.
[Obrazek: https://lx-net.pl/walet/image/motylek1.jpg]
Algorytm ten zwany jest również motylkowym lub radix-2 ze względu na graficzne podobieństwo zobrazowania par obliczeń do skrzydeł motyla. Na każdym kroku algorytmu łączy się do obliczeń ze sobą pary wartości stąd wymaganie aby ilość danych wejściowych była potęgą liczby 2.
A teraz nareszcie w skrócie na temat samego algorytmu complex fft zastosowanego w kodzie Waleta.
[Obrazek: https://lx-net.pl/walet/image/fftc_dsp.jpg]
To co poprzednio oznakowałem jako blok cfft w rzeczywistości składa się z dwóch głównych części. Pierwsza, skojarzona z podzielonymi przez MultiOutput danymi I oraz Q zwana iqReceiver z klasy CallbackStream pozwala na odczyt próbki 2*256 danych z kanałów I oraz Q a drugą część to właściwa complex fft składająca się z dwóch procedur: dsps_fft2r_fc32 oraz dsps_bit_rev_fc32. Na końcu uzyskuje się macierze wartości rzeczywistych i urojonych dla każdego bin-a, które służą do wyliczenia jego amplitudy na zasadzie składania wektorów: globalFullAmplitudes[k] = sqrtf(real_part * real_part + imag_part * imag_part). MultiOutput zapewnia dostarczenie danych I oraz Q oraz uruchomienie funkcji obliczeń transformaty po skompletowaniu 2*256 danych poprzez procedurę writeDataCallback.
W pliku tekstowym wydzieliłem najważniejsze polecenia dotyczące complex fft
Analiza danych z tego pliku pewnie zainteresuje tylko programistów ale muszę przyznać, że uzyskanie tego kodu jest na granicy moich możliwości programistycznych, szczególnie przy obsłudze danych i procedur dotyczących dużych struktur oraz funkcji i metod Arduino Audio Tools. Posiłkowałem się przy tym AI, która oczywiście halucynowała proponując nieraz dziwne rozwiązania i przepraszając za błędy ale w końcu udało się uzyskać poprawny wynik.
Ostateczna procedura complex fft nie jest częścią Arduino Audio Tools ale stanowi natywne opracowanie producenta procesora stąd na początek ładowana jest biblioteka #include "esp_dsp.h" a w dalszej części uruchamiane jest sprzętowe wspomaganie fft.
W wyniku działania opisanych procedur powstaje tabela globalFullAmplitudes zawierająca 256 binów ułożonych według następującego porządku: bin 0 to tzw. składowa stała, która w tym przypadku jest pomijana, bin 1 to wartość dla 93.5Hz, bin 2 to wartość dla 187Hz ..... bin 127 to wartość dla 12000Hz, bin 128 to wartość dla -12kHz, ...... bin 254 to wartość dla -187.5Hz i bin 255 to wartość dla częstotliwości -93.75Hz. Jak widać dane w tabeli są dla częstotliwości dodatnich i ujemnych podane w odwróconym porządku co należy uwzględnić w graficznym obrazowaniu tych danych.
Udało mi się również ostatnio uzyskać obrazowanie sygnału z pomocą fft również w formie "wodospadu" ale o tym napiszę w kolejnym poście jak przygotuję zarówno rysunki, zdjęcia oraz nagrania i obiecuję, że nie będzie więcej tak mocno programistycznych postów, choć oczywiście postaram się coś wyjaśnić w razie pytań i potrzeby.
Szykuję też pliki Arduino i binarne dla Github ale również potrzebuję trochę czasu na uporządkowanie kodu przez jego prezentacją.
cdn. lj
Github