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

symulacyjny model filtru polifazwoego

cepa 03 Lip 2004 11:37 67704 164
  • #31
    h-doc
    Poziom 27  
    shg - filtry IIR mogą być alternatywą, ale jest jeden problem - tutaj potrzebna jest większa precyzja obliczeń. O ile FFT możesz policzyć stałoprzecinkowo na 16-bitach, o tyle filtr IIR bałbym się robić z taką precyzją. Niestety IIR są bardzo wrażliwe na precyzję obliczeń i nie dość, że charakterystyka może ulegać fluktuacjom, to w najgorszym wypadku filtr się wysypie. Dlatego IIR powinno się robić na zmiennym przecinku, a to w przypadku 8-bitowych stałoprzcinkowych procków mogłoby być trudne.
    Jest jeszcze inny pomysł (zastosowany w kodeku MPEG). Mianowicie filtry polifazowe. Wymagają zarówno sporo pamięci jak i mocy obliczeniowe, ale są stabilne (bo są na bazie FIR), a co również istotne - umożliwiają odtworzenie sygnału pierwotnego (rzecz niemożliwa przy DFT czy FFT).
  • IGE-XAO
  • #32
    shg
    Specjalista techniki cyfrowej
    Tak, racja zapomniałem :oops: W tym PDFie z załącznika są wzory do obliczenia potrzebnej liczby bitów (tam jest to chyba liczone na stałym przecinku, ale głowy nie dam sobie uciąć). Im niższa częstotliwość tym więcej bitów trzeba, tam dla najniższych wystarczą chyba 24bity (o ile mnie pamięć nie myli, bawiłem się tym dawno temu ostatecznie przedłem na FIR).
  • #33
    hunterhouse
    Poziom 26  
    jednak jak bym chciał zrobić to w bascomie to i tak i tak muszę liczyć na zmiennym przecinku bo bascom funkcje trygonometryczne liczy tylko na zmiennych typu single
    więc sprubuje zrobić to przez DFT i IIR i zobacze co działa szybciej

    a tak w ogule to mam jeszcze pytanie czy dobrze rozumiem to IIR
    najpierw muszę zaprojektować analogowy filtr wedłóg schematu 1 w pdf kolegi

    następnie muszę wyliczyć alfa , beta i gama które są stałe dla każdego filtru

    a następnie muszę postępować według algorytmu z rysunku 9a

    ??????????????????????????????????????????????????? [/img]
  • #34
    shg
    Specjalista techniki cyfrowej
    hunterhouse napisał:
    najpierw muszę zaprojektować analogowy filtr wedłóg schematu 1 w pdf kolegi
    następnie muszę wyliczyć alfa , beta i gama które są stałe dla każdego filtru
    a następnie muszę postępować według algorytmu z rysunku 9a

    Dokładnie tak :!: , ale musi być jakieś ale :D . Podane zostały jedynie uproszczone wzory na alfa, beta i gamma, wykorzystujące SAA (Small Angle Aproximation). To taka zasada, że dla małych kątów sin(x) jest w przybliżeniu równe x, wobec czego zakłada się, że sin(x)=x (to zostało tam gdzieś napisane i przez to radykalnie upraszczają się wzory, podobnie założono tan(x)=x). Dla częstotliwości f<fs/8 (fs-częstotliwość próbkowania) wzory te są jak najbardziej poprawne i w miarę dokładne. Jednak przy wyższych częstotliwościach zaczynają się schody. Np. przy fs=44.1 kHz, projektując filtry na 12 i 16 kHz (przy użyciu SAA) okazuje się, że mają prawie takie samo pasmo przepustowe :D Sprawdzałem to na MathCADzie. Wzorów na wyższe pasma nie podali, trzeba je samemu wyprowadzić, ale to też nie wiele pomaga, bo wyprowadziłem je (straszne potworki), a parametry poprawiły się tylko nieznacznie. Przy wysokich częstotliwościach, nie tylko trzeba skorzystać z pełnych wzorów, ale jeszcze wziąć poprawkę na przesunięcie fazowe i kwantyzację w domenie czasu, a otym jak to zrobić już nic nie piszą, podają tylko wyniki (Sprawdziłem. Wyniki moich obliczeń bez uproszczeń odbiegały od tych podanych w tabeli dla 8 i 16 kHz). Dlatego olałem IIR i zająłem się FIRem, ale FIR w BASCOMIE? jakoś mi się nie widzi, IIR jest na pewno szybszy. Na zmiennym przecinku, analizator widma, w bascomie? odświeżanie raz na 10 sekund gwarantowane. Wielką zaletą IIRów jest to, że po obliczeniu współczynników (możesz obliczyć je na kompie i wrzycić do programu jako tabelę, bo jak już wiesz są stałe) nie trzeba wykonywać skomplikowanych obliczeń (tylko dodawanie i mnożenie). DFT natomiast wymaga albo ogromnej ilości funkcji trygonometrycznych, albo tabeli. Robiłem kiedyś na Amidze (a to dużo szybsze niż taki Atmel mały) DFT w asemblerze na stałym przecinku. Odświeżanie - masakra, jakieś 5 razy na sekundę.

    Rozważania nad sposobem realizacji analizatora:

    Jeżeli ten analizator widma nie ma być przyrządem pomiarowym to z powodzeniem wystarczy 16bitowa precyzja, (chyba 48 cykli na mnożenie i 4 na dodawanie/odejmowanie dwóch liczb 16bitowych).
    Przy prędkości zegara=24MHz, cykl trwa 0.5us.
    na jedną próbkę - 3 mnożenia (144 cykle), 4 dodawania/odejmowania (16cykli) + inne pierdoły (ze 20 cykli)=180 cykli (po drobnej optymalizacji) = 90us/próbkę
    90us/próbkę*10 kanałów=900us/próbkę.
    dodajmy ze 100 us na wyświetlanie itp.
    wychodzi 1ms na próbkę.
    czyli częstotliwość próbkowania 1kHz.
    Wniosek:
    tak się nie da zrobić analizatora widma, chyba że na szybszym procu.

    Inne podejście:
    założenia j.w.
    Będziemy próbkować do RAMu i dopiero przetwarzać.
    Minimalna częstotliwość na analizatorze = 30Hz Maksymalna - 10kHz, potrzeba próbkowania min. 20 kHz
    żeby zmieścić jeden okres sygnału 30 Hz potrzeba 20000/30=666 :D próbek
    czas akwizycji: 1/30 s
    1 ms/ próbkę = 666 ms na cały bufor, czyli to już by mogło działać, ale odświeżanie 1.5 razy na sekundę to trochę mało. no itrzeba mieć odpowiednią ilość pamięci, a w małych Atmelach z tym ciężko.

    Inne podejście:
    Podział pasma na wysokie i niskie
    założenie jest takie, żeby przetwarzać równolegle dwa strumienie A i B
    strumień A o wysokiej czestotliwości próbkowania, ale w małym buforze (chodzi o rozmiar w czasie :D )
    strumień B o niskiej częstotliwości próbkowania w dużym buforze (również czas)
    fizycznie oba bufory mogą mieć ten sam rozmiar (w bajtach tym razem)
    Strumień B można by uzyskać przez decymację (wycinanie co którejś próbki) ze strumienia A. Pozostaje jednak niekorzystne zjawisko aliasingu. Jednakże sygnały o małych częstotliwościach mają zazwyczaj znacznie większą amplitudę, od sygnałów o dużych częstotliwościach, więc aliasing przy większym poziomie sygnału byłby niewielki.
    można też równolegle (stworzyć 3 strumień) odfiltrować wyższe składowe strumienia A (tylko czym IIRem ? :D ) i dopiero dokonać decymacji. To by działało, tylko trzeba znaleźć jakiś szybki filtr dolnoprzepustowy (FIR i DFT odpadają). Właściwie, to można dorzucić jedenastego IIRa, skoro już i tak jest ich 10 :D
    No to dzielimy góra (A) - 1kHz i wyżej (do 10kHz), dół (B) - mniej niż 1kHz (do 30Hz)
    A i C - próbkowanie 20 kHz
    bufor na A - 20000/1000=20 bajtów :D
    bufor na C (strumień równoległy, tyle, że odfiltrowany) - 20 bajtów
    strumień B próbkowanie 1 kHz
    Bufor na B - 1000/30=33 bajty
    Właśnie wymyśliłem szybki filtr dolnoprzepustowy - wystarczy uśrednić strumień A - stała czasowa 1ms, akurat pasuje :D 20 dodawań 8bit+16bit jakieś 70 cykli. dzielenie 16bit przez 20 (20 - stała, więc istnieje możliwość jakiejś optymalizacji) - nie wim ile, zakładam 50 cykli. razem 120 cykli = 60 us Tak uśrednioną próbkę odjąć od wszystkich próbek w strumieniu A, żeby wyeliminować aliasing (tylko, że taki odwrotny), czyji wyciąć składową stałą. 20 cykli (10us), razem 70us
    70 us+ 5 IIRów na strumieniu A * 90 us/próbkę * 20 próbek = 9070 us na cały strumień A. Załóżmy, że 10ms
    co 10 ms dodawana próbka do bufora B, zapełniony w 330ms
    5 firów * 90us/próbkę * 33 próbki = 14ms na obliczenia
    Całkowity czas całkowitego odświeżania :D :
    330+14+10=354ms, czyli prawie 3 razy na sekundę. Ale wysokie pasmo (A) można odświeżać znacznie częściej, bo aż 100 razy na sekundę :D. można by jeszcze zastosować na wejściu audio kilka analogowych filtrów , multiplekser, i przetwornik A/C, ale to już nie jest prawdziwe DSP :D ....
    A co by było gdyby antyaliasing (filtrację) zrobić sprzętowo? Np. tylko 2 filtry, jeden dolno-, drugi górnoprzepustowy, multiplekser sterowany z proca, częste próbkowanie do pasma A (sygnał analogowy przez filtr górnoprzepustowy), i od czasu do czasu (powiedzmy po zapełnieniu bufora A) przełączał by się na pasmo B (sygnał analogowy przez filtr dolnoprzepustowy) i pobierał jedną próbkę do bufora B? Tego dowiecie się w następnym pasjonującym odcinku "DSP i Transmisja", a tym czasem kalkulator do filtrów IIR:
    http://www.dsptutor.freeuk.com/IIRFilterDesign/IIRFilterDesign.html

    Miłych rozważań nad filozoficznymi aspektami Nieskończonej Odpowiedzi na Impuls, życzy kierownictwo ośrodka w Tworkach. :D
  • IGE-XAO
  • #35
    shg
    Specjalista techniki cyfrowej
    :arrow: h-doc (post z 23 czerwiec 2004):
    przedstawiony przez kol. hunterhouse algorytm DFT w BASICU nie oblicza połowy, tylko całe widmo, i przedstawia je w postaci rzeczywistej (składowych sin i cos)przy wykorzystaniu IDFT można na podstawie tego widma odtworzyć sygnał (więc na 100% widmo jest kompletne). dwa, to piszesz, że po filtracji da się odtworzyć sygnał, a po FFT się nie da (wracając do tego kodeka).... ja nie rozumiem, bo przecież po FFT można odtworzyć bez problemu, przy użyciu IFFT.
    tu masz link:
    http://www.xilinx.com/ipcenter/catalog/logicore/docs/vfft32.pdf
    strona 3 wzór 3 - FFT, wzór 4 - IFFT.

    polecam też: www.dspguide.com, chyba najlepsza książka o DSP, w dodatku za darmo!

    :arrow: hunterhouse
    Wielkość prążka widma po DFT oblicza sie tak:
    wszystki elementy tablic re[] i im[] dzieli się przez N, a potem dla każdego
    sqtr(re[x]*re[x]+im[x]*im[x]), ale po angielsku to sie nie nazywa amplitude, tylko magnitude (czyli wielkość(?)).
  • #36
    shg
    Specjalista techniki cyfrowej
    Opowieści o IIR część druga, niekoniecznie ostatnia.

    System wygląda następująco: (od wejścia do wyjścia)
    - wejście audio
    - filtry dolno, i górnoprzepustowy na 1 kHz
    - multiplekser analogowy
    - przetwornik A/C
    - procesor (sterujący też multiplekserem)
    - weświetlacz

    od strony programowej:
    dwa strumienie A i B, jak w części pierwszej.
    strumień A - próbkowanie 20 kHz / 20 bajtów bufora
    strumień B - próbkowanie 1 kHz / 33 bajty bufora
    Analiza:
    Przełączamy multiplekser na strumień A (przez filtr górnoprzepustowy),
    pobieramy 20 próbek z częstotliwością 20 kHz (czas: 1ms)
    Przełączamy multiplekser na strumień B (przez filtr dolnoprzepustowy), pobieramy 33 próbki z częstotliwością 1 kHz (czas: 33 ms)
    "Odpalamy" IIRy na strumieniu A
    90us/próbkę * 20 próbek * 5 kanałów = 9ms
    "odpalamy" IIRy na strumieniu B
    90us/próbkę * 33 próbki * 5 kanałów = 14.8ms
    Suma czasu przetwarzania:
    1ms + 33ms + 9ms + 14.8 ms = 57.8 ms na kompletne widmo, co daje odświeżanie 17 razy na sekundę
    Jeżeli coś pokaszaniłem, to wybaczcie, ale pierwszego posta pisałem w środku nocy (a nawet rano, dla niektórych), a drugiego oparłem na założeniach pierwszego. dziś przeczytałem wszystko jeszcze raz dokładnie, i wydaje mi się sensowne :D. Miłej lektury.

    P.S. Chyba sam sobie zrobie taki analizator :D
  • #37
    morph13
    Poziom 25  
    Przyznam, że opis jest "łopatologiczny", ale znakomicie wyjaśniający dany temat (zjawisko), nareszcie konkretni ludzie na forum, to dobrze wróży elektrodzie w przyszłości........ :-). Sam przyznam, że opis wyjaśnił mi kilka niewiadomych :oops: .
  • #38
    shg
    Specjalista techniki cyfrowej
    Wiadomość z ostatniej chwili :!:
    Nie można uruchamiać IIR na takich krótkich próbkach, bo na początku i na końcu ramki IIR jest niestabilny, trzeba odrobinę wydłużuć bufory.
    Następnie zrobić coś, co się nazywa fade in i fade out, czyli stopniowe zwiększanie amplitudy sygnału na początku bufora i stopniowe zmniejszanie na końcu. Wszystko wygląda tak samo jak wcześniej pisałem, ale po próbkowaniu jest jeszcze tyle mnożeń (na stałym przecinku 8bit ofkoz), ile próbek w buforze. Nie wiem jak bardzo należy wydłużyć bufory, myślę że razy 2 powinno bez problemów wystarczyć, no ale im krótszy bufor, tym czas przetwarzania krótszy. Myślę, że odświeżanie nie powinno spaść poniżej 8-10 razy na sekundę.
    A te filtry wejściowe powinny mieć strome charakterystyki, można zastosować filtry Czebyszewa nawet o dość sporych zafalowaniach w pasmie przepustowym, ludzkie zmysły są i tak dość ułomne :D .
  • #39
    h-doc
    Poziom 27  
    Riposta:
    1. a jednak ten programik w Basicu liczy tylko połowę widma - porównaj pętlę z linii 90 (liczymy od zera do N/2) ze wzorem na DFT (liczymy od zera do N-1)
    2. odnośnie FFT i IFFT - wzory to niestety nie wszystko. To, że sygnał poddany transformacie FFT a następnie IFFT (bez żadnych dodatkowych modyfikacji) wróci do postaci pierwotnej jest oczywiste. Jednak kto robi FFT bez okienkowania i zakładkowania? Uwzględnij jeszcze przeciek widma oraz błedy powstające na skutek skończonej precyzji obliczeń a odpowiedź sama się napisze... Poza tym jeśli byłoby to takie proste, to wskaż mi jakiś kodek, który koduje widmo a następnie na jego podstawie potrafi odtworzyć sygnał oryginalny.
    Powtarzam - nie jest możliwe odzyskanie sygnału pierwotnego z jego widma (chyba, że tak jak pisałem - robimy to bez nakładek, okienek i liczymy z dużą precyzją)
  • #40
    shg
    Specjalista techniki cyfrowej
    :arrow: h-doc. A czy kolega potrafi odróżnić transformatę na liczbach zespolonych, od transformaty na liczbach rzeczywistych?
    Wzór podany przez kolegę jest na liczbach zespolonych, na wyjściu otrzymujemy tyle samo liczb co na wejściu.
    Kod w BASICU jest na liczbach rzeczywistych i na wyjściu też otrzymujemy tyle samo liczb co na wejściu. Co z tego, że liczy do N/2 ?, wstawia widmo do dwóch tablic, i liczb jest tyle samo.
    Teoria Nyquista (w skrócie): Pasmo przenoszenia jest 2 razy mniejsze, niż częstotliwość próbkowania. Przykład. próbkowanie 2kHz, czas 1s, mamy 2000 próbek. zgodnie z teorią Nyquista mamy pasmo 1 kHz. Zapuszczamy program kol. hunterhouse (N=2000) i otrzymujemy widmo w zakresie 0-N/2Hz z rozdzielczością 1Hz.
    Jeżeli dalej uważasz, że to połowa widma to pokaż mi skąd mogę ściągnąć (albo zrób sam) sygnał o częstotliwości 1.5kHz próbkowany z częstotliwością 2kHz :D

    >Jednak kto robi FFT bez okienkowania i zakładkowania?
    Ja robię :D Napisałem analizator widma i nie stosuję okienkowania, bo to tylko zniekształca widmo, a ja muszę widzieć je dokładnie takie, jakie jest (zakładkowania jednak używam).

    Jakiś kodek, który odtwarza sygnał z widma - JPEG
    transformata DCT, też w pewien sposób tworzy widmo sygnału :D, potem koder poddaje je kwantyzacji, i kompresji bezstratnej
    dekoder dekompresuje i przepuszcza przez IDCT (odwrotna dyskretna transformata cosinusów) i tak oto z widma odtwarza się oryginalny sygnał :D . Kilka nieistotnych szczegółów procesu (de)kodowania pominąłem.
    Do odtworzenia sygnały z widma wcale nie jest potrzebna super precyzja. Jeżeli sygnał składał się z ośmiobitowych próbek, to można obliczyć jego widmo, zapisać je w postaci również ośmiobitowych próbek i odtworzyć na jego podstawie sygnał oryginalny bez strat. (sygnał i widmo zajmują tyle samo miejsca). Proponuję dokształcić się w tej dziedzinie, bo to co napitałem, to są podstawy podstaw :D
  • #41
    h-doc
    Poziom 27  
    shg - po pierwsze nie wymądrzaj się (to odnośnie Twojej propozycji mojego dokształcania);
    po drugie - kolego, liczba punktów transformaty FFT jest zawsze równa liczbie próbek, które poddano transformacie (i to są chyba te "podstawy podstaw", o których raczyłeś wspomnieć). Zatem jeśli próbek masz N, a transformatę liczysz do N/2 to jak widać liczysz tylko połowę widma. Kolega wspominając coś o liczbach rzeczywistych i zespolonych coś kompletnie pokręcił. Gwoli wyjaśnienia dla mniej zorientowanych - widmo sygnału rzeczywistego wykazuje symetrię względem połowy środka zakresu. Stąd też wystarczające dla sygnałów rzeczywistych jest liczenie połowy widma (do N/2). Nie zmienia to jednak faktu, że jest to połowa widma i można tak robić tylko dla sygnałów rzeczywistych. Przy liczeniu FFT twierdzenie Nyquista obchodzi mnie tyle co zeszłoroczny śnieg - przetwarzając próbki nie muszę wiedzieć kompletnie nic o częstotliwości próbkowania (gdzie w algorytmie FFT masz częstotliwość próbkowania???)
    po trzecie - no cóż, skoro robisz FFT bez okienkowania, to nieźle wiarygodne musisz uzyskiwać wyniki :-)
    po czwarte - DCT to nie FFT. Poza tym jak już wspomniałeś JPEG to kompresja stratna (i to nie tylko na skutek kwantyzacji - również na skutek transformacji DCT/IDCT). Oprócz tego z sygnałem wizyjnym postępuje się zupełnie inaczej niż z sygnałem audio (w JPEG transformacji DCT poddaje się próbki bez nakładkowania i okienkowania, więc sygnał nie jest tak zniszczony, jak w przypadku FFT na sygnale audio).
    po piąte - w dalszym ciągu czekam na (rozsądny) algorytm, który najpierw zakoduje sygnał przy pomocy FFT, a potem odtworzy go przy pomocy IFFT, mogący działać na bieżąco w sposób ciągły (czyli bez dostępu na raz do wszystkich próbek).
  • #42
    shg
    Specjalista techniki cyfrowej
    h-doc napisał:
    shg - po pierwsze nie wymądrzaj się (to odnośnie Twojej propozycji mojego dokształcania);

    Dobra, przepraszam, troche przesadziłem. :oops: . Nie mam nic na swoją obronę, oprócz tego że byłem niewyspany (jak zwykle) i ciężkiej książki, o której niżej. Ten post nie jest żadną docinką z mojej strony, bronię tylko swoich poglądów i wyjaśniam kilka spraw popierając wszystko materiałem dowodowym najprzedniejszego sortu :D . A nóż (widelec) się komuś przyda.

    h-doc napisał:
    Kolega wspominając coś o liczbach rzeczywistych i zespolonych coś kompletnie pokręcił.

    Co pokręciłem ?
    Obie (o tym niżej) transformaty opierają się na tym samym wzorze na liczbach zespolonych, ale:
    Wzór podany przez kolegę jest dla zespolonych danych wejściowych i wyjściowych.
    Program w BASICU to numeryczna implementacja tego algorytmu (czyli wszystko liczone na rzeczywistych, bo na kompie inaczej się nie da), ale okrojona do liczb rzeczywistych na wejściu.
    Otworzyłem sobie The Scientist and Engineer's Guide to Digital Signal Processing" rozdział 8 strona 147. (w wersji elektronicznej :D )
    http://www.spectrumsdi.com/ch8.pdf (476kB)
    Wyraźnie tam napisano, że są dwie DFT - na liczbach rzeczywistych i zespolonych. I co jescze? Do zespolonych wzór taki, jak podał kolega, do rzeczywistych, taki jak w programie w BASICu.
    Otworzyłem też rozdział 31:
    http://www.innerlighttheory.com/ch31.pdf (389kB, na innym serwerze, dziwne....)
    I przeczytałem, co następuje:
    DSPguide napisał:

    Each of these four transforms has a complex version and a real version. The
    complex versions have a complex time domain signal and a complex frequency
    domain signal. The real versions have a real time domain signal and two real
    frequency domain signals. Both positive and negative frequencies are used in
    the complex cases, while only positive frequencies are used for the real
    transforms.
    oraz:
    if the frequency domain is real and periodic (rzezcywista, okresowa, czyli jakaś ramka z danymi, to jest w PDFie dokładnie opisane przyp. by ja), the
    appropriate section is over the one-half cycle containing the positive
    frequencies, either 0 to N (continuous frequency domain (ta nas nie obchodzi przyp. by ja)) or 0 to N/2 (discrete
    frequency domain (ta natomiast nas dotyczy przyp. by ja)).

    Czyli wychodzi na to, że dla liczb rzeczywistych (a na takich operujemy) odrzuca się UJEMNE częstotliwości, bo są niepotrzebne i zostaje tylko ta "połówka", która jest kompletnym widmem dla sygnałów rzeczywistych
    Przeglądałem kilka implementacji FFT i w wielu widziałem, że były osobne procedury dla liczb zespolonych i rzeczywistych te dla rzeczywistych dawały na wyjściu 2 tablice o rozmiarz N/2, a te dla zespolonych, jedną tablicę liczb zespolonych o rozmiarze N.

    h-doc napisał:
    widmo sygnału rzeczywistego wykazuje symetrię względem połowy środka zakresu. Stąd też wystarczające dla sygnałów rzeczywistych jest liczenie połowy widma (do N/2).

    Widmo będzie miało lustrzane odbicie, ale tylko wtedy, gdy jako dane dla DFT na liczbach zespolonych poda się dane z zerową składową urojoną, czyli liczby rzeczywiste, bo jak już napisałem DFT dla rzeczywistych liczy tylko do N/2.

    h-doc napisał:
    Nie zmienia to jednak faktu, że jest to połowa widma i można tak robić tylko dla sygnałów rzeczywistych.

    A skąd wziąć zespolony sygnał audio? Z przetwornika A/C (i ogólnie w cyfrowych systemach audio) jest tylko rzeczywisty. (no dobra są jakieś tam analizatory stosowane w energetyce itp., ale tu chodzi o przetwarzanie dżwięku).

    h-doc napisał:
    Przy liczeniu FFT twierdzenie Nyquista obchodzi mnie tyle co zeszłoroczny śnieg - przetwarzając próbki nie muszę wiedzieć kompletnie nic o częstotliwości próbkowania (gdzie w algorytmie FFT masz częstotliwość próbkowania???)

    Mnie częstotliwość próbkowania nic nie obchodzi, to był tylko przykład, równie dobrze mogło być wszystko razy 29 :D .
    Twierdzeniem Nyquista też się nie przejmuję, bo zawsze liczę do N/2 :D .
    Chodziło mi o to że dla sygnału rzeczywistego jest o połowę mniej prążków w widmie, niż próbek sygnału, w związku z czym nie da się spróbkować sygnału o okresie krótszym niż 2 próbki, a tak jest właśnie w przypadku 1.5kHz próbkowanego z fs 2kHz
    Chyba teraz mi się udało przetworzyć myśli na text :D )

    h-doc napisał:
    po trzecie - no cóż, skoro robisz FFT bez okienkowania, to nieźle wiarygodne musisz uzyskiwać wyniki :-)

    Całkiem wiarygodne!
    Tu masz screena z mojego analizatora widma (13kB)
    http://shogoonn.w.interia.pl/elektronika/gfft1.gif
    Dwie częstotliwości - 1kHZ + 1.3kHz, próbkowanie 28kHz 8bit, bez okienkowania i bez zakładkowania, transformata 1024 punkty (FFT na liczbach rzeczywistych).
    Te "zbocza" (ang. spectral leakage) to efekt aliasingu (rzecz całkowicie normalna, bo częstotliwości sygnału nie trafiają dokładnie w częstotliwości analizatora), postrzępiony przebieg przy wyższych częstotliwosciach, to szum kwantyzacji (przy 8 bitach nawet go słychać :D).
    Tylko z takiego "nieobrobionego" widma można odtworzyć sygnał.
    Po najechaniu myszą na fragment widma analizator podaje częstotliwość, nutę (!) i poziom sygnału (to jest screen ze starszej wersji, bez poziomu sygnału). Już niejednego nauczyłem się grać na gitarze dzięki niemu :)

    h-doc napisał:
    JPEG to kompresja stratna (i to nie tylko na skutek kwantyzacji - również na skutek transformacji DCT/IDCT).

    Też mi stratność, jeżeli ustwisz kompresję na 0 (nie mówię tu o bezstratnym formacie JPEG, bo to zupełnie co innego), to DCT/IDCT bierze i tak udział w procesie (de)kodowania a straty są nie do zauważenia.
    Sprawdziłem to tak (poziom strat, nie to, czy dane przeszły przez DCT/IDCT, bo tego jestem pewien po lekturze źródeł kodeków JPEG):
    wygenerowałem obrazek (szum, żeby miał dużą entropię, niestety nie biały, bo takowego w Corelu nie znalazłem), zapisałem jako BMP i jako JPEG z kompresją ustawioną na 0 (JPEG był większy :).
    tu masz linki do obrazków:
    http://shogoonn.w.interia.pl/elektronika/jpeg.jpg [100k]
    http://shogoonn.w.interia.pl/elektronika/bmp.bmp [76k]
    http://shogoonn.w.interia.pl/elektronika/delta.bmp [76k]
    W Corelu 9 (przy okazji znalazłem BUGa - przy wklejaniu obiektu, nie wstawia go w pozycję 0,0, tylko trochę obok, musiałem ustawić ręcznie, z tym też był problem :D ) wykonałem na nich operację różnica (efekt w pliku delta.bmp)
    Znakomita większość (71968) pikseli nie różni się od oryginału, natomiast reszta (4832) różni się o 1 bit (LSb ofkoz). Strata prawie żadna, no a dla obrazu, czy dżwięku nie ma to kompletnie żadnego znaczenia.
    Przy okazji - niezła metoda na sprawdzenie stratności przy kompresji :D.
    Co innego text, programy, no ale kompresji transformacyjnej się tu nie stosuje, bo i tak by nic nie dała.
    Mając na myśli odtwarzanie sygnału z jego widma chodziło mi właśnie o coś takiego, czyli coś, czego na oko/ucho odróżnić się nie da od oryginału.
    A tak jeszcze do tej precyzji, to np. JPEG nie zapisuje widma na zmiennym przecinku, tylko jako bajt (podczas przetwarzania, zaraz po DCT), a po kwantezacji zostaje z tego kilka bitów :D . Są wprawdzie implementacje zmiennoprzecinkowe, ale nie wiem po jaką cholerę, bo i tak kwantyzacja (ta w procesie kompresji) daje więcej strat, niż zmienny przecinek zysku :).

    h-doc napisał:
    (w JPEG transformacji DCT poddaje się próbki bez nakładkowania i okienkowania, więc sygnał nie jest tak zniszczony, jak w przypadku FFT na sygnale audio).

    Echhhh. A kto ci każe robić zakładkowanie i okienkowanie? Te dwa procesy stosuje się jedynie do wizualizacji danych, aby usunąć zbocza powstałe w wyniku aliasingu (okienkowanie) i rozciągnąć widmo w czasie, co poprawia jego czytelność (zakładkowanie). A że niektórzy tępo powielają je wszędzie gdzie się da, nawet tam gdzie nie powinno się tego robić, to już inna historia.
    Znalazłem kiedyś equalizer na FFT zdaje się na sourceforge.net
    Działał tak: wejście->FFT->mnożenie przez odpowiednie wagi (zgodnie z ustawieniami suwaków)->IFFT->wyjście
    Bez zakładek, bez okienek. (z racji powolności FFT wymagał procka 800MHz żeby chodził w czasie rzeczywistym, zdaje się :D)

    Tutaj masz krótkie info o okienkowaniu, jak widać służy ono wyłączne do wyeliminowania "spectral leakage" pzy wizualizacji sygnału (strona baaaaardzo znanego producenta sprzętu pomiarowego :D ):
    http://www.tek.com/Measurement/cgi-bin/framed..._Notes/fft/window.html&FrameSet=oscilloscopes

    h-doc napisał:
    po piąte - w dalszym ciągu czekam na (rozsądny) algorytm, który najpierw zakoduje sygnał przy pomocy FFT, a potem odtworzy go przy pomocy IFFT, mogący działać na bieżąco w sposób ciągły (czyli bez dostępu na raz do wszystkich próbek).

    A z tym to już nie do mnie, ja się zajmuję stricte przetwarzaniem, a nie kompresją (niby też przetwarzanie :D), więc za bardzo na kodekach się nie znam (oprócz JPEG, którego usiłowałem zaimplementować w asemblerze na Amidze, nic z tego nie wyszło, ale takie były moje początki :D ).
    Chyba, że mam napisać w basicu program, który policzy widmo (rzeczywiste), a następnie z widma odtworzy oryginalny sygnał, bo to się da zrobić :D . Tylko powiedz mi o co chodzi z tym działaniem w sposób ciągły, bo bufor na próbki i tak trzeba mieć.

    Post sprzed miesiąca:
    lesiotech napisał:
    Chodzi tu, podobnie jak w przypadku formatu MPEG-Video, o najefektywniejszą w chwili obecnej metodę kompresji.

    Z tego co mi wiadomo, to obecnie najefektywniejsza metoda kompresji wykorzystuje transformatę falkową (np. popularny pod linuchami vorbis).
    Vorbis 48kbps STEREO brzmi lepiej niż MP3 64kbps MONO, chociaż przy tym bitrate to i tak shitowo :D to moje subiektywne odczucie.

    Posta zmieniałem parę razy, bo w czasie pisania czytałem The Scientist and Engineer's Guide to Digital Signal Processing" i nie tylko, w zasadzie potwierdziło się to co wiedziałem i trochę dokształciłem się, jeżeli chodzi o zespoloną transformatę :D.

    Przepraszam, że ten post taki długi, ale inaczej pisać nie umiem, przynajmniej się człowiek nie nudzi, bo jest co poczytać :D .

    Nie chcę wojny na posty, bo chyba nie po to jest elektroda. (dobra wiem, że sam zacząłem :( ).

    Proponuję kontynuować (tak, tu jest napisane kontynuować, nie zamknąć :D ) temat, bo jest naprawdę ciekawy.

    Ze swej strony mam pytanie odnośnie filtrów IIR:
    Czy ktoś wie jak obliczać współczynniki, przy f bliskim fs/2 :?:
  • #43
    h-doc
    Poziom 27  
    No, w końcu doszliśmy do konsensusu odnośnie DFT. Standardowy sygnał audio jest rzeczywisty i można liczyć połowę widma. Należy jednak o tym pamiętać, że jest to tylko połowa widma (bo można się nieźle na tym przejechać jeśli zacznie się operować na liczbach zespolonych).

    Zespolony sygnał audio? Proszę bardzo - zajmujesz się przecież przetwarzaniem - pewnie więc robiłeś transformatę Hilberta albo modulację kwadraturową - tu już trzeba przedstawiać sygnał w postaci zespolonej.

    Aha, tak odnośnie FFT zrobionej tylko do N/2 (dla sygnału rzeczywistego). Zrób sobie potem IFFT z tego okrojonego sygnału. Ciekawe, co otrzymasz...

    Odnośnie próbkowania sygnału 1.5kHz z częstotliwością 2kHz to nie widzę problemu - z pomocą przychodzi podpróbkowanie.

    odnośnie http://shogoonn.w.interia.pl/elektronika/gfft1.gif - Ty to nazywasz wiarygodnym widmem? Dla mnie nie ma tu dwóch prążków, tylko bardziej przypomina to szum o odpowiednio ukształtowanej charakterystyce, a powinny być 2 wyraźne prążki i cisza w okolicy. Proszę Cię - zrób sobie to samo dla kilku wybranych okien czasowych, a wyniki zdecydowanie się poprawią. A tak na marginesie - pod czym pisałeś ten soft - Amiga?

    Te przykłady z JPEG były dosyć interesujące - o to chodzi w JPEG, że oko ludzkie nie jest specjalnie wrażliwe na tego typu zniekształcenia. Poza tym weź pod uwagę, że JPEG nadaje się do obrazów naturalnych. Z obrazami syntetycznymi (np. tekst, wykresy spod Excela) jest już dużo gorzej. Ciekawe, jak wypadłyby Twoje testy dla szumu barwnego.

    Odnośnie zmiennoprzecinkowego przetwarzania w JPEG - cóż, są to na tyle złożone algorytmy, że łatwiej jest tutaj operować na zmiennym przecinku. Podobnie jest zresztą z MPEG - norma też jest zmiennoprzecinkowa, chociaż oczywiście istnieją stałoprzecinkowe implementacje. Jest to jednak trudniejsze do zrobienia, bo trzeba borykać się z ograniczeniem zakresu, nasyceniami itd.

    Twoje wątpliwości odnośnie okienkowania i zakładkowania są ciekawe, ale chyba jednak te operacje są potrzebne. Znowu wracam do MPEG. W kodeku audio MPEG wykorzystuje się FFT w modelu psychoakustycznym (gdzie nie ma przecież wizualizacji danych, tylko przetwarzanie). I niestety - stosowane tam jest dosyć głębokie zakładkowanie i okienkowanie sygnału, a rezultaty są przecież całkiem dobre.

    Odnośnie projektowania filtrów IIR - z czym masz dokładnie problem? Czy metoda prototypu analogowego nie zdaje egazaminu?

    pozdrawiam!
  • #44
    shg
    Specjalista techniki cyfrowej
    Taaa soft piszę na Amidze :D . (dawniej 600, obecnie 1200)
    Te nieostre prążki wcale mi nie przeszkadzają, bo widzę i tak jaka to częstotliwość, a tylko to mnie interesuje :D, a założenie było takie: Dawno, dawno temu, nie maiałem dostępu do netu, więc nie mogłem sobie ściągać tabów na gitarę, pomyślałem, że trzeba jakoś wykombinować nutki samemu :D Ze słuchu grać umiem, ale mam ograniczenie prędkości percepcji, więc napisałem ten analizator. Uruchamiam analizator, wrzucam do niego mój ulubiony kawałek, kilkadziesiąt minut ślęczę nad widmem i już wiem jak to zagrać :D . Ale tu masz rację, to jest wizualizacja widma i warto zastosować okienkowanie. Przy solówkach to widmo jeszcze jako tako wyglądało, ale jak grał cały zespół, to słabsze dźwięki tonęły w "wyciekach widma". Pisałem program w czasach, gdy dopiero zaczynałem zabawę z DSP i nie przejmowałem się okienkami, a do tego żeby zrobić zakładkowanie doszedłem sam :D .

    Zmiennoprzecinkowość JPEG - oczywiście standard został zdefiniowany i wszystkie obliczenia były na zmiennym przecinku, ale zaraz znaleźli się tacy, którzy chcieli szybciej i więcej i zaimplemantowali to na stałym przecinku. Na Amidze np. jest tylko jeden (słownie: 1) dekoder JPEGa zmiennoprzecinkowy, reszta na stałym i to z mnóstwem optymalizacji. Różnica, aż w oczy kole: zmiennoprzecinkowy (na Ami bez koprocesora numerycznego) jest około 100 razy wolniejszy od najszybszego obecnie FastJPEG :D

    Filtry IIR
    Wszystko działa ładnie i pięknie, ale tylko dla niskich częstotliwości, w momencie gdy częstotliwość jest wyższa od fs/5, filtry zachowują się dziwnie, tak jak pisałem wcześniej, ich pasmo przepustowe przestaje zgadzać się z pasmem wynikającym z projektu (nie ma mowy o małej precyzji, bo wszystko liczę na typie double (MathCAD)), konkretniej, to częstotliwość środkowa tak obliczonego filtru jest niższa od założonej :( .
    W tym PDFie, co wrzuciłem są wzory, na ich podstawie (musiałem przekształcać, ale o tym pisałem wcześniej - problem z SAA) wyliczyłem
    współczynniki i kapa, dla niskich f, owszem działa, ale dla wysokich środek pasma filtru wypada niżej (dużo niżej).
    Prawdopodobnie problemem są przesunięcia fazowe i nieuniknione straty procesu kwantyzacji w domenie czasu. Motorola (bo to ich PDF) podaje tam tabelę z gotowymi współczynnikami, te z kolei zachowują się całkiem przyzwoicie i filtry mają pożądaną charakterystykę, tylko nie napisali szuje jak je policzyli, A ja muszę liczyć wszystko od nowa, bo potrzebuję różnych częstotliwości próbkowania, a oni wyliczyli tylko dla 44.1kHz.
    Przejrzałem kilka projektów equalizerów na IIR, wszystkie mają już wyliczone współczynniki (w tabeli) dla 44.1kHz, niektóre też dla kilku innych częstotliwości, ale żaden nie potrafi wyliczać ich dla dowolnej częstotliwości próbkowania, ani innych parametrów filtru. :(
    Czytałem, szukałem, sprawdzałem różne sposoby obliczeń, ale wszystkie dają podobne wyniki - marne :(
    Jeżeli coś nie jasno opisałem, to mogę wrzucić screeny z MathCADa, tam to ładnie widać.
    Owszem istnieje rozwiązanie tego problemu - oversampling, ale jak się ma procesor 14MHz (w architekturze x86 odpowiednik - 386 28MHz, może nawet trochę szybciej), to nie ma mowy o takim rozwiązaniu.
  • #45
    hunterhouse
    Poziom 26  
    zrobiłem próbnie na początek samo serce analizatora według algorytmu IIR i czciałbym się zapytać czy nie ma tam żadnych błędów


    Dim A(10) As Single 'tablica alfa
    Dim B(10) As Single 'tablica beta
    Dim C(10) As Single 'tablica gama
    Dim X(3) As Byte 'tablica próbek z a/c
    Dim Y(3) As Byte 'tablica wyników w obliczeniach
    Dim I As Byte 'I
    Dim J As Byte 'J
    Dim K As Byte 'K

    Dim Y1(10) As Byte 'tablica y1 dla różnych częstotliwości
    Dim Y2(10) As Byte 'tablica y1 dla różnych częstotliwości
    Dim Y3(10) As Byte 'tablica y1 dla różnych częstotliwości
    Dim Wyn(10) As Byte 'tablica wyników do wyświetlenia
    Dim Wyn_p As Word 'zmienna pomocnicza przy obliczeniacz
    Dim Tmp As Single 'zmienna pomocnicza przy obliczeniacz
    Dim Tmp_2 As Single 'zmienna pomocnicza przy obliczeniacz
    Dim Nr As Byte 'określa który filtr liczymy


    Config Timer1 = Timer , Prescale = 1
    'dane do tablic alfa , beta i gama
    A(1) = 0.000787462865
    A(2) = 0.00157244917
    A(3) = 0.00316016172
    A(4) = 0.00628062774
    A(5) = 0.0124054279
    A(6) = 0.0242101804
    A(7) = 0.461841095
    A(8) = 0.0845577687
    A(9) = 0.1199464
    A(10) = 0.159603

    B(1) = 0.498425074
    B(2) = 0.496855102
    B(3) = 0.493679677
    B(4) = 0.487438745
    B(5) = 0.475189144
    B(6) = 0.451579639
    B(7) = 0.407631781
    B(8) = 0.330884463
    B(9) = 0.2601072
    B(10) = 0.1800994

    C(1) = 0.998415336
    C(2) = 0.996816209
    C(3) = 0.993522095
    C(4) = 0.986812425
    C(5) = 0.972715729
    C(6) = 0.941937749
    C(7) = 0.871031797
    C(8) = 0.699565951
    C(9) = 0.3176087
    C(10) = -0.4435172

    Main:
    'pętla główna
    Do
    Gosub Sample

    For Nr = 1 To 10
    Gosub Iir
    'out wyn(nr) - wyświetlenie na diodach
    Next Nr
    Loop



    Iir:
    Timer1 = 0
    Start Timer1
    'oblicza IIR dla jednej częstotliwości zależnej od "nr".

    For I = 1 To 3

    'pobiera interesujące nas wyniki z poprzednich obliczeń.
    Y(1) = Y1(nr)
    Y(2) = Y2(nr)
    Y(3) = Y3(nr)

    'ustala wartośc "i", "j", "k".
    Select Case I
    Case 1
    J = 2 : K = 3
    Case 2
    J = 3 : K = 1
    Case 3
    J = 1 : K = 2
    End Select

    'oblicza według wzoru.
    Tmp = X(i) - X(j)
    Tmp = Tmp * A(nr)
    Tmp_2 = C(nr) * Y(k)
    Tmp = Tmp + Tmp_2
    Tmp_2 = B(nr) * Y(j)
    Tmp = Tmp - tmp_2
    Tmp = Tmp + Tmp

    'zapisuje wynik
    Select Case I
    Case 1 : Y1(nr) = Tmp
    Case 2 : Y2(nr) = Tmp
    Case 3 : Y3(nr) = Tmp
    End Select
    Next I

    Wyn_p = Y1(nr) + Y2(nr)
    Wyn_p = Wyn_p + Y3(nr)
    Wyn(nr) = Wyn_p / 3

    Stop Timer1
    Cls
    Lcd Timer1 ; "clk "
    Wait 2
    Return




    Sample:
    'próbkuje 3 próbki z częstotliwością 44.1 khz
    X(1) = 124
    X(2) = 230
    X(3) = 54
    Return

    dodałem też timer żeby sprawdzić ile zajmuje liczeni i wychodzi jakieś od 7000 do 11000 cykli zegarowych a to i tak dobrze bo zegar jest 16Mhz
    używam uładu atmega8515 chociaż nie ma przetwornika bo taki akurat mam w domu i nie mogę przetestować w wersji real
  • #46
    h-doc
    Poziom 27  
    shg - no nieźle - nie wpadłbym na to, że można zastosować analizator FFT do rozgryzania solówek gitarowych - podziwiam Twoją pomysłowość!
    Jeśli chodzi o filtry IIR, mogę Ci podrzucić metodę bazującą na prototypie analogowym. Najpierw musisz zaprojektować normalny filtr analogowy. A potem możesz pójść dwiema drogami - albo transformacja biliniowa albo metoda niezmienności odpowiedzi impulsowej. To dosyć dobre metody projektowania filtrów IIR. Jeśli jednak chciałbyś pozotać przy swojej metodzie, koniecznie sprawdź, czy bieguny transmitancji zawierają się w okręgu jednostkowym (warunek stabilności).
  • #47
    shg
    Specjalista techniki cyfrowej
    :arrow: h - doc
    Z prototypem analogowym jest ten ból, że mój program będzie musiał projektować dodatkowo filtr analogowy.
    Chcę zrobić to metodą opisaną w PDFie. Przeczytałem go jeszcze raz, tam są wszystkie potrzebne równania !!!
    Mi to nie działało, bo nie uwzględniłem czegoś przy wyprowadzaniu wzorów, w tej chwili nie pamiętam czego, bo to było dawno. Oni nawet napisali jak to policzyć, ale widoczne wtedy nie zauważyłem, gotowego rozwiązania jednak nie dali, trzeba wyprowadzać wzory-potwory. Dziś już próbowałem, MathCAD meczył się pół minuty i wywalił coś szerokiego na 2 ekrany (czcionka 9, 1280x960). Spróbuję jeszcze raz wyprowadzić, na pewno da się prościej.
    Jeżeli to, co masz to kod źródłowy, to bardzo chętnie przejrzę, bo samej matematyki znalazłem więcej niż dużo :D , ale nie chce mi się przez nią "przegryzać". Zresztą w tym PDFie jest właśnie prototyp analogowy, a potem wszystko tak sprytnie poprzekształcane, żeby podstawić f środkową i dobroć.
    I jeszcze jedno, najlepiej, jakby metoda projektowania dawała 3 współczynniki (alfa, beta gamma), bo wszystko, co widziałem daje 6 (filtr pasmowoprzepustowy), albo 12 (filtry dolno- i górnoprzepustowy połączone w kaskadę). Algorytm na 3 współczynnikach jest o wiele szybszy od tego na 6, czy 12, a w tym wypadku szybkość odgrywa kluczową rolę.
  • #48
    shg
    Specjalista techniki cyfrowej
    :arrow: hunterhouse
    Dobrze, ale nie wszystko, nie zgadza się z algorytmem :(
    W sumie to się zgadza, ale tylko dla 3 próbek, ale te 3 próbki to za mało.

    aby zapisać sygnał 31 Hz (najniższa częstotliwość filtrowana) przy próbkowaniu 44.1kHz potrzebujesz aż: 44100/31=1423 próbki. Do atmega można chyba podłączyć pamięć zewnętrzną (ale nie wiem, bo nigdy nie używałem atmegi), wtedy by działało, ale pozostaje i tak problem szybkości :( .
    Przetworzenie takiej ilości danych zajmie prockowi duuużo czasu. Szczególnie, że obliczenia wykonujesz na zmiennym przecinku. W tym wypadku konieczne jest stosowanie stałego przecinka.
    Jak by procek był na tyle szybki, że przetwarzał by dane w czasie rzeczywistym, to niepotrzebny był by bufor (oprócz tych 3 próbek, wymaganych przez algorytm.)
    pętla liczącza IIR , będzie działać tylko dla 3 próbek, a mamy więcej, więc:
    Code:
    i = 1 ' początkowe i
    
    for n = 1 to 3
     x(n)=0 : y(n)=0
    next
    wyn(nr)=0 ' dlaczego tak to pozniej
    For n = 1 To liczba_probek ' "n", a nie "i", bo i liczy 1,2,3,1,2,3....
     x(i)=probka ' w tym momencie zostaje pobrana kolejna probka z bufora (pozycja n w buforze), lub przetwornika A/C (jeżeli działa w czasie rzeczywistym)
    'ustala wartośc "i", "j", "k".  - nie "i", bo to już mamy ustalone
     ...
    'oblicza według wzoru.
     ...
     tmp=y(i)
    ' i tu pojawia się problem, bo dane wyjściowe filtru zajmują tyle samo, co dane wejściowe, a że filtrów jest 10, potrzeba by na nie ponad 14kB pamięci.
    ' należy zostawić to, co jest i na bieżąco liczyć amplitudę na wyjściu, wtedy bufory wyjściowe nie będą potrzebne. W ŻADNYM WYPADKU NIE WOLNO UŚREDNIAĆ SYGNAŁU: średnia dąży do 0, jeżeli sygnał nie zawiera składowej stałej, a po przepuszczeniu przez te filtry nie zawiera. Najprościej będzie uśredniać wartość bezwzględną, dla analizatora widma wynik będzie zadowalający.
     wyn_p=abs(tmp) ' tmp - to co zostało zapisane do y(i) - czyli wyjściowy sygnał filtru
     wyn(nr)=wyn(nr)+wyn_p ; sumowanie wartosci bezwzglednych
     i=i+1
     if i>3 THEN i = 1
    next
     wyn(nr)=wyn(nr)/liczba_probek ' suma/ilość = średnia

    Ideałem było by obliczenie wartości skutecznej (RMS) sygnału, co nie jest takie skomplikowane :D
    Wartość skuteczna to pierwiastek kwadratowy ze średniej podniesionych do kwadratu wartości próbek :D

    więc zamiast tego:
    Code:

     wyn_p=abs(tmp) ' tmp - to co zostało zapisane do y(i)
     wyn(nr)=wyn(nr)+wyn_p
     i=i+1
     if i>3 THEN i = 1
    next
     wyn(nr)=wyn(nr)/liczba_probek ' suma/ilość = średnia

    Trzeba dać to:
    Code:

     wyn_p=tmp*tmp ' podnosi do kwadratu
     wyn(nr)=wyn(nr)+wyn_p ' sumuje
     i=i+1
     if i>3 THEN i = 1
    next
     wyn(nr)=wyn(nr)/liczba_probek ' suma/ilość = średnia
     wyn(nr)=sqrt(wyn(nr)) ' pierwiastek i mamy RMS

    RMS spowolni obliczenia niezauważalnie, bo i tak są już powolne, wszystko przez zmienny przecinek.
    Jak zauważełeś pozbyłem się tablic wyjściowych dla poszczególnych filtrów (Y1() itd.), wszystko jest w tablicy y(), bo i tak na raz liczony jest tylko jeden filtr i zaraz po tym otrzymujemy gotowy do wyświetlenia wynik.
    Pisałem o niestabilności IIR na początku i końcu ramki. Tu nie należy przejmować się tym problemem, bo ramka jest dość duża w stosunku do rozmiaru filtra i ewentualne zniekształcenia spowodowane niestabilnością będą niezauważalne.
    Postaram się jak najszybciej zaimplementować to na stałym przecinku (jak tylko rozpracuję metodę obliczania filtrów), można spróbować zrobić to na 16 bitach, przy tej rozdzielczości minimalna częstotliwość to nieco ponad 50Hz, zgodnie z tym co piszą w PDFie, ale jeszcze trzeba by się upewnić, w którym miejscu postawić przecinek, bo pisali coś że wartości mogą przekroczyć 1, wtedy filtr się wyłoży :(
    Wciąż pozostaje problem szybkości, proponuję jednak filtrację analogową i podział na 2 strumienie danych, tak jak pisałem w jednym z moich wywodów, przyspieszenie jak widać radykalne :D . Ja przynajmniej tak zrobię.
    Analizator będę robił na 89Cx051 24MHz, większość kodu będzie w asemblerze, ale postaram się najpierw całość zrobić w BASICu, żeby było bardziej zrozumiale, jak coś będzie w asm, to postaram się dobrze skomentować :D
    A i jeszcze bardzo ważna rzecz - dane wejściowe to liczby zmiennoprzecinkowe z zakresu (-1.0 do 1.0) może być i większy zakres, ale lepiej zostawić taki.
    Nie wiem, czy wszystko, co tu napisałem jest dobrze, bo nie sprawdzałem w praktyce, ale wygląda poprawnie :D .

    Jeszcze jedno, tak mi teraz do głowy wpadło:
    Czy ktoś z osób śledzących ten wątek (szczególnie koledzy hunterhouse i h-doc) zna następujące języki programowania: C, PASCAL i Modula2 ?
    Pytam, bo napisał bym najpierw oprogramowanie w języku AmigaE, który jest bardzo podobny do wyżej wymienionych i żródła udostępniał bym na bieżąco, szybciej były by wyniki :D .
    Odpowiedzi proponuję na prywatną wiadomość, żeby nie zaśmiecać forum.
  • #49
    hunterhouse
    Poziom 26  
    napisałem też algorytm DFT, ale w bascomie działa dużo wolniej
    jak się postaram to może wyciągnę 10 x na sekundę ale procesorem z 16Mhz zegarem

    dzieli pasmo na dwie częstotliwości i próbkuje osobno i liczy osobno
    w sumie 2x32 próbki
    próbkuje z f=32kHz i 2kHz

    muszę kupiś sobie proca z a\c i troche diodek, wtedy sprawdzić jak to wychodzi w praktyce
    jak się zdecyduję na DFT to się poświęce i napiszę w asm to będzie dużo szybsze
    natomiast IIR i zmienno-przecinkowców to w asemblerze nie widze

    tak dla informacji oprócz bascoma umiem pascala i asm(oczywiście)
  • #50
    shg
    Specjalista techniki cyfrowej
    DFT z podziałem pasma, to też niezły pomysł, na pewno będzie łatwiej, bo nie trzeba liczyć współczynnkiów. No i przede wszystkim, łatwo to zrobić na stałym przecinku.
    Wyświetlanie - raczej koniecznie w skali logarytmicznej, w liniowej będzie kaszanka, bo już robiłem na Amidze :D
    Prosty, szybki i nie wymagający tablic sposób na skalę logarytmiczną:
    Wynik do wyświetlenia masz w postaci 16 bitowej bez znaku, wyświetlacz 16 LEDów.
    Sprawdzasz po kolei od MSb czy bit jest 1, czy 0,
    jak 0, to następny,
    jak 1, to zapalasz 16 minus tyle LEDów ile było bitów na 0

    W asm to tylko kilka instrukcji:
    Code:

    MOV R0, #0 ' na razie góra nie świeci
    MOV A, MSB ' najbardziej znaczacy bajt do acc
    JZ policz_dol ' jeżeli zero, to wyświetla tylko dolne ledy
    MOV R0, #8 ' licznik pozycji
    petla_msb:
     RLC  A   ' wstawia kolejne bity z acc do C, zaczynając od MSb
     JC policz_dol
     DJNZ R0, petla_msb

    policz_dol:
    MOV R1, #0 ' dół nie świeci
    MOV A, LSB ' najmniej znaczacy bajt do acc
    JZ wysw
    MOV R1, #8 ' licznik pozycji
    petla_lsb:
     RLC  A
     JC wysw
     DJNZ R1, petla_lsb

    wysw:
     ... 'zapala tyle ledów z górnej ósemki, ile jest w R0
     ... 'zapala tyle ledów z dolnej ósemki, ile jest w R1
     RET

    Pisałem "z głowy", ale powinno być dobrze, da się to pewnie jeszcze zoptymalizować, pod konkretny sposób wyświetlania
  • #51
    hunterhouse
    Poziom 26  
    mam jeszcze takie pytanie czy w DFT do analizatora widma konieczne jest liczenie tablicy z wynikami imx (to tablica zawierająca część urojoną)

    narazie wygląda to tak
    jeden plik to wersja pierwsza
    a dróga (-asm) przyspieszona ale opiera się na tech samych założeniach ma wstawki asemblerowe i skasowane niepotrzebne działania
  • #52
    shg
    Specjalista techniki cyfrowej
    Niestety trzeba obie tablice policzyć, chyba że.... (o tym niżej) .
    Pomimo tego, że wejściowe dane są rzeczywiste, to na wyjściu są już zespolone i zapełniają obie tabilie, więc Im[] nie jest pusta. Oczywiście można by spróbować olać jedną z tablic Re[], albo Im[], obojętnie.
    Mogło by to doprowadzić do dziwnego zachowania analizatora, jeżeli spełnione zostały by odpowiednie warunki fazowe, pomimo obecności sygnału na wejściu, analizator mógł by niczego nie wyświetlać. Przy sygnale audio prawdopodobieństwo małe, więc można i warto spróbować.
  • #53
    h-doc
    Poziom 27  
    shg - programuję w C (no i czasem we wszelkiej maści asemblerach)
    Jeśli chodzi zaś o projektowanie filtrów IIR przy pomocy prototypu analogowego, to mam pełne algorytmy - do policzenia ręcznego lub w MathCad. W C też by to można policzyć (chociaż rozbijanie na ułamki proste jest nieco kłopotliwe). Oczywiście z prototypu dolnoprzepustowego można przejść na filtr górno- lub pasmowo-przepustowy.
    Jeśli chodzi o samo wykonanie filtru, to jest już kwestia jak potraktujesz ułamek przedstawiający transmitancję (może to być pojedynczy filtr - jeden ułamek albo na przykład grupa sekcji bikwadratowych - suma ułamków stopnia drugiego).
  • #54
    shg
    Specjalista techniki cyfrowej
    :arrow: hunterhouse
    trochę błędów było, przejrzałem tylko wersję w BASICU.
    Nie tak jak trzeba liczyłęś kąty (chyba, nie wiem, od razu poprawiłem, tak żeby było jak najszybciej i oszczędniej).
    Poprawiłem podprogram liczący sin() i cos().
    Tabela sin(), nie musiała być aż taka wielka. Pełen okres powinien zawierać się w tylu próbkach, ile masz próbek sygnału. Tutaj tabela Sin() ma 32 wartości dla pełnego okresu. Można zrobić tak jak próbowałeś, tabelę z jednej ćwiartki (8 wartości), ale sama procedura licząca sin i cos zajmie więcej miejsca, niż uzyskasz przez skrócenie tabeli, poza tym jak jest cała tabela, to jest szybciej. Można pójść o krok dalej i na końcu tabeli powtórzyć pierwszą ćwiartkę, wtedy będzie można wywalić jedno AND 31 (to przed liczeniem Tmp_c, po dodawaniu (tmp_x=tmp_x+8).
    Wyeliminowałem możliwość powstania błędów OVERFLOW.

    Tablica dane zawiera 32 próbki ośmiobitowe ze znakiem (-128 do 127)
    Widmo wynikowe (wielkości prążków) są w tablicy Rex() - liczby z zakresu
    0-256 - stałoprzecinkowe, można było uzyskać większą dokładność, ale trzeba by trochę pokombinować.
    Dzielenia należy zastąpić przesuwaniem bitów w prawo, problem jest taki, że nie wiem, jak to się ma do BASCOMA, czy robi przesunięcie logiczne, czy arytmetyczne (ze znakiem). Coś mi się jednak zdaje, że logiczne, nie wiem bo nie sprawdzałem.

    Program skompilowałem, olałem to, że tablice "wyjechały" poza pamięć. obliczenie całej transformaty zajmuje 420ms na 89C2051 24MHz
    Kod jest baaardzo nieoptymalny, ale przez to chyba przejrzysty. :D

    Wpliku jest tylko procedura DFT, funkcje liczące sin i cos i tabela
  • #55
    shg
    Specjalista techniki cyfrowej
    h-doc - mogą być pliki z MathCADa, tylko czy to mi wygeneruje współczynniki takie, że bedzie działało poprawnie dla wysokich częstotliwości? Znalazłem coś, miało działać tak jak trzeba, ale działało tylko lepiej :( . Wciąż nie mogę uzyskać odpowiednich charakterystyk dla 8 i 16 kHz przy fs 44.1kHz :(

    MAM!
    Ściągnąłem DSP extensiopn pack do MathCADa, są tam funkcje liczące współczynniki IIR, Filtry działają tak jak powinny.
    Tylko co ja mam z tymi współczynnikami zrobić?
    Nie znalazłem nic o topologii filtru, doszedłem tylko do tego, że to filtry dolno- i górnoprzepustowy połączone kaskadowo.
    Niema też nic o jego równaniu :(
    Jedyne czego się dowiedziałem, to że pierwsza kolumna macierzy jest licznikiem, a druga mianownikiem.

    Może te zrzuty ekranu z MathCADa coś komuś podpowiedzą:
    Tu jest obrazek z parametrami i wyliczonymi współczynnikami, sygnał do przetworzenia (biały szum), jego widmo i wzmocnienie filtru w zależności od f
    http://shogoonn.w.interia.pl/elektronika/img/mcd1.png [6kB]
    Widmo przefiltrowanego szumu i dla porównania jeszcze raz charakterystyka filtru:
    http://shogoonn.w.interia.pl/elektronika/img/mcd2.png [4kB]

    Wszystkie algorytmy, które są w DSP Extension Pack przedstawiono w książce Signal processing algorithms autorzy - Samuel Stearns i Ruth David, ale co mi po tym, nidostępna w wersji elektronicznej :(

    :!: Zmieniłem program do wyświetlania w skali logarytmicznej, znalazłem błąd, poprawiłem, a przy okazji samo się trochę zoptymalizowało :D :!:

    MathCAD wylicza takie same współczynniki, jak ten kalkulator, do którego podałem link :D Czyli można powiedzieć, że już znalazłem :D .
  • #56
    h-doc
    Poziom 27  
    Faktycznie, wygląda na to, że pierwsza kolumna to współczynniki licznika transmitancji a druga to współczynniki mianownika.
    Żeby uzyskać charakterystykę filtru, lepiej podać impuls Diraca niż szum biały.
  • #57
    shg
    Specjalista techniki cyfrowej
    :arrow: h-doc
    Dobra wiem, już jak wygląda równanie transmitancji filtru :D.

    A równanie dla programistów jest takie:
    Code:
    y[n]=a0*x[n]+a1*x[n-1]+a2*x[n-2]-b1*y[n-1]-b2*y[n-2]
    

    a0 do a2 to pierwsza kolumna, b0-b2 to druga.
    kilka uproszczeń:
    a1 jest zawsze równe 0,
    a0=-a2
    mamy więc:
    y[n]=a0*(x[n]-x[n-2])-b1*y[n-1]-b2*y[n-2]
    i jewszcze jedna ciekawa rzecz. Pamiętacie program w BASICU z PDFa?:
    y[n]=2*(alfa*(x[n]-x[n-2])+gamma*y[n-1]-beta*y[n-2])
    Uderzające podobieństwo :D

    Kalkulator do IIR, nieco lepszy(?) niż poprzedni, generuje kod źródłowy w C i nie wykłada się przy wyższych f (sprawdzałem):
    http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html

    W pierwszym kalkulatorze IIR (kilka postów wcześniej), jest kod źródłowy w Javie, wylicza poprawne współczynniki, więc teoretycznie wszystko, co potrzeba już mamy.

    Jeszcze raz screenshoty z MathCADa, tym razem impuls Diraca, ale filtr już nieco inny.

    http://shogoonn.w.interia.pl/elektronika/img/mcd1.png [7kB]
    http://shogoonn.w.interia.pl/elektronika/img/mcd2.png [5kB]

    Od razu ładniej wygląda. :D
    Zabiorę się za robotę jak tylko zrobię preemfazę do nadajnika :D .
  • #58
    hunterhouse
    Poziom 26  
    mam jeszcze pytanie dotyczące danych z przetwornika w jakiej postaci muszą one byś zapisane (ile bitów i czy ze znakiem)


    ten twój program dla całego widma avr morze zrobić 13 razy na sekundę ale dochodzi jeszcze 2x próbkowanie i wyświetanie multiplekserem w międzyczasie
  • #59
    shg
    Specjalista techniki cyfrowej
    8 bit powinno wystarczyć, koniecznie ze znakiem.
    Większość przetworników daje niestety sygnał bez znaku, jeżeli sygnał będzie 8 bit bez znaku, to wystarczy odjąć od niego 128 i już jest ze znakiem.
    Przy ośmiobitowym sygnale poziom szumu w widmie to ok -72dB, czyli do analizatora zpokojnie wystarczy.

    :arrow: h-doc
    w kodeku MPEG audio jest chyba jakaś FFT z "rozwiniętymi" pętlami?
    Przeglądałem kiedyś źródła pod linuxem i gdzieś coś takiego widziałem.
    Może dało by się to wrzucić zamiast DFT, szybkość była by piorunująca :D . Z tego co pamiętam ta transformata była 32, albo 64 punktowa, postaram się jeszcz przejrzeć to, co mam, może znajdę.
    Tak baj de łej, ściągnąłem sobie jakieś FFT, na Athlonie 2100+ XP przetworzenie 5 min. muzyki z CD (transformata 1024 punkty) zajmuje mniej niż sekundę :D .
  • #60
    h-doc
    Poziom 27  
    no, teraz ten filtr wygląda tak jak trzeba.
    Jeśli chodzi o FFT w MPEG to z tego co pamiętam to w modelu psychoakutycznym 1 wykorzystywana jest 512-punktowa transformata a w modelu 2 transformata jest 1024-punktowa. Jeśli chodzi o samo FFT, to w źródłach jest już pewna optymalizacja algorytmu, ale nie wgryzałem się w nią specjalnie - znacznie większym obciążeniem była filtracja polifazowa - i tutaj trzeba było przysiąść, żeby w miarę szybko chodziła.