logo elektroda
logo elektroda
X
logo elektroda
Adblock/uBlockOrigin/AdGuard mogą powodować znikanie niektórych postów z powodu nowej reguły.

symulacyjny model filtru polifazwoego

cepa 05 Sty 2003 11:47 69522 164
  • #1 76162
    cepa
    Poziom 12  
    Znalazłem proste i dobrze działające funkcje fft i ifft. Za ich pomocą robie analize 512 próbek. Na wyjściu mam 2 tablice liczb typu double. Spektrum uzyskuje w ten sposób (wyciete ze znalezionego kodu):

    double *real, ///czesc rzeczywista
    *imag, ///ta druga
    *spec; ///widmo

    ///...

    fft(512,real,imag);
    double a = 0;
    int i = 0;
    for(i = 0; i < 512; i++)
    {
    spec[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);
    if(spec[i] > a) a = spec[i];
    }
    for(i = 0; i < 512; i++)
    spec[i] /= a;

    ///...

    I teraz mam problem. Wszystko da sie wyswietlic tylko ze zamiast poukładanego spektrum, mam sieczke, tzn: zadna liczba nie jest w tym miejscu w którym powinna byc. Gdy zwieksze czestotliwosc sygnału podanego w real to widmo sie przesówa ale nie pokolei tylko skacze z jednego konca na drugi zblizając sie ku srodkowi.

    Wiec wie ktos jak z czegos takiego uzyskac ładne poukładane widmo do wyswietlenia na wykresie ??
  • #2 430042
    morph13
    Poziom 25  
    Prawdopodobnie dane wejsciowe powinny być poddane funkcji zmieniajacej położenie bitów w kazdej próbce. Ksiązkowo opisywane to jest jako "motylek". Jak się dokop[ie do literatury to coś podeślę.
  • #3 432081
    h-doc
    Poziom 27  
    trzeba dokonać tzw. inwersji bitowej, która uporządkuje prążki.
  • #5 675302
    lesiotech
    Poziom 13  
    możesz coś napisać ifft?
    projektujr radiostacje cyfrową z wykorzystaniem kompresji dżwięku,Kompresja ma poleać na analizie widma,i wycince najsłabszych częstotliwości.
    do dekompresji potrzebuje ifft a nie znalazłem do tej pory alorytmu:(
    Czy jak na wejście algorytmu fft podam spektrum to na wyjściu uzyskam sygnał??
  • #6 676126
    h-doc
    Poziom 27  
    lesiotech - bardzo mi przykro, ale metoda kompresji w dziedzinie częstotliwości chyba nie zda egazaminu. Trzeba kompresować w dziedzinie czasu. Można to zrobić tak jak w mpeg, czyli przy pomocy filtrów polifazowych. A może wogóle zaimplementować mpega (np. I warstwę, która jest najprostsza)?
  • #7 679085
    marmur99
    Poziom 17  
    Z tymi bitami to jest tak, że je przestawiasz, np:

    000 -> 000
    001 -> 100
    010 -> 010
    011 -> 110
    100 -> 001
    itd.

    czyli próbka z pozycji 1 teraz będzie na pozycji 4 a próbka z pozycji 2 pozostaje na miejscu.

    Swego czasu pisałem to na kompie w C++ i moja procedurka porządkująca wygląda tak:

    void Porzadkuj()
    {
    int dlugosc_pozycji = log10(ilosc_probek)/log10(2);
    int nowa_pozycja;
    complex tmp;

    for (int i=1;i<ilosc_probek;i++)
    {
    nowa_pozycja=0;
    for (int o=0;o<dlugosc_pozycji;o++)
    {
    int bit = (i & (int(pow(2,(dlugosc_pozycji-1-o)))));
    if (bit!=0) nowa_pozycja+=int(pow(2,o));
    }
    if (nowa_pozycja>i) //zapobieganie dwukrotnej zamianie miejsc
    {
    tmp=tab_dane[i];
    tab_dane[i]=tab_dane[nowa_pozycja];
    tab_dane[nowa_pozycja]=tmp;
    }
    }
    }

    Oczywiście na mikrosterownik będzie to wyglądało nieco inaczej, ale ogólna idea jest taka.

    Pozdrawiam,

    Marmur99
  • #8 679738
    h-doc
    Poziom 27  
    marmur99 - ja z kolei robiłem to metodą tabelaryczną, co było dużo szybsze (choć zajmowało nieco pamięci).
  • #9 679935
    lesiotech
    Poziom 13  
    h-doc napisał:
    lesiotech - bardzo mi przykro, ale metoda kompresji w dziedzinie częstotliwości chyba nie zda egazaminu. Trzeba kompresować w dziedzinie czasu. Można to zrobić tak jak w mpeg, czyli przy pomocy filtrów polifazowych. A może wogóle zaimplementować mpega (np. I warstwę, która jest najprostsza)?

    Format MPEG-Audio (Wydanie II) zbudowany jest podobnie jak algorytm PASC (Precision Adapptive Subband Coding) firmy Philips, stosowany na kasetach DCC (Digital Compact Cassete - cyfrowa kaseta kompaktowa). Chodzi tu, podobnie jak w przypadku formatu MPEG-Video, o najefektywniejszą w chwili obecnej metodę kompresji. Punktem wyjścia dla omawianej metody jest model psychoakustyczny. Jak we wszystkich algorytmach kompresji pracujących ze stratami, wykorzystuje się tu niedostatki możliwości percepcyjnych człowieka. Z przetwarzanego sygnału usuwa się nie tylko informacje redundantne (nadmiarowe), lecz także częstotliwości które nie będą zauważone. Przy użyciu 32 filtrów pasmowych analizuje się spektrum częstotliwości akustycznych z podziałem na subpasma. Sygnały leżące poza zakresem słyszalności człowieka nie są uwzględniane. Wykorzystuje się przy tym efekt maskowania - sygnały głośniejsze zagłuszają sąsiadujące z nimi sygnały słabsze, tak więc te słabsze są maskowane i nie muszą być rejestrowane. Może się wiec zmieniać długość słowa bitowego.

    sorryże nie na temat,le wydaje mi sie że podstawą mpega jest właśnie analiza częstotliwości,i jeszcze jakieś bzdety,które dla prostoty chce pominąć bo mi sie Atmega ie wyrobi.
    [/u]
  • #10 680092
    h-doc
    Poziom 27  
    Brawo kolego lesiotech! Tak się składa, że na pracę dyplomową robiłem kodek MPEG2 audio. W MPEG-u jest analiza widma FFT. Jest ona wykorzystywana właściwie we wszystkich modelach psychoakustycznych. Natomist sygnał, który przechodzi przez koder (nie przez model psychoakustyczny) nie podlega analizie FFT, tylko (jak sam podkreśliłeś) filtracji podpasmowej. Jest to (jak sama nazwa wskazuje) filtracja, a nie transformata. W związku z czym sygnał można po takiej filtracji odtworzyć (w przeciwieństwie do sygnału po analizie FFT).
    Poza tym w przytoczonym przez Ciebie cytacie jest drobny błąd (zresztą mocno rozpowszechniony) - chodzi o "nie rejestrowanie" sygnałów maskowanych. Otóż po raz kolejny powtarzam - kodek MPEG nie "wycina" maskowanych dźwięków tylko je kwantyzuje. Tylko w skrajnych przypadkach (np. przy małej przepływności bitowej) podpasmo zostaje całkowicie wyzerowane.

    Chcę jeszcze dodać, że algorytm MPEG audio (nawet najprostsza I warstwa) jest to rozbudowany i czasochłonny wynalazek. Mój kodek chodził ledwo w czasie rzeczywistym, przy mocno wyżyłowanej optymalizacji (wstawki w asemblerze, większość rzeczy liczona równolegle i na stałym przecinku). Robiłem to na procku Infineon Carmel DSP z zegarem 20MHz.

    Pozdro!
  • #11 680333
    morph13
    Poziom 25  
    Może nie wniosę tym postem nic sensownego, ale ten temat idealnie pasuje do "DSP i Transmisja". Może admin się tym zajmie :-)

    A z mojej strony jeszcze dodam, że fajnie gdy pojawiają się wstawki z procedurami, jest co porównywać no i temat się rozwija.
  • #12 681689
    lesiotech
    Poziom 13  
    Dzięki kolego h-doc z acene informacje,Czy w takim razie mogę prosić o jakieś sugestie??
    Zależy mi na kompresji dźwięku(mowy) do 8 kbit/s próbkowanie 8KHz i 8 bitów co bez kompresji daje 64kbit/s
    wydaje mi sie że mpeg to w tym przypadku lekka przesada ;)

    Przyznaje że nie mam w tej kwestii zbyt dużej wiedzy ale za to algorytmizacja i pisanie w asm to dla mnie bułka z masłem. :D
  • #13 681751
    h-doc
    Poziom 27  
    pytanie zasadnicze - jakie zasoby sprzętowe masz do dyspozycji? Do MPEGa musi być proc sygnałowy albo pecet, bo zwykły mikrokontroler nie ma wystarczającej mocy obliczeniowej.
    Są jeszcze kodeki gsm, ale one sprawdzają się jedynie przy sygnale mowy (kodowanie muzyki raczej odpada).
    Może wystarczy kompresja DPCM albo ADPCM? Myślę, że to by się dało zrobić na mikrokontrolerze, chociaż nie wiem, czy 8-bitowa precyzja wystarczy. Przy 8 kHz próbkowania, i przy zegarze 16 MHz masz do dyspozycji 2000 cykli na każdą próbkę, więc jest to całkiem sporo.
    Spróbuj sobie coś skompresować przy użyciu ADPCM i zobacz, czy rezultaty są zadowalające, a potem pomyślimy nad ewentualną implementacją.
  • #14 681791
    lesiotech
    Poziom 13  
    do dysppozycji mam atmege 168 -20Mhz i 16kB flash
    może być 16 bitów ale efekt końcowy to bezwzględnie nie więcej niż 8kbit/s

    jak uda mi się dobyć at91m55800 to bede miał 100MHz 16 bitowego procka(porównywalny z 486SX 100Mhz) wiec mpega hyba pociągnie
    zobacze z adpcm i pomyśle o inych rozwiązaniach.
  • #15 695341
    piotr_go
    Konstruktor DIY elektronika
    patrzcie co ostatnio przypadkiem znalazłem
    http://instruct1.cit.cornell.edu/courses/ee476/FinalProjects/s1999/lisa/
    "Spectrum analyzer" na 8535
    może komuś sie przyda :)
  • #17 697287
    hunterhouse
    Poziom 26  
    mam takie pytanie
    chciałem zrobić analizator widma oparty na FFT ale niestety nieumiem C
    więc pytam się czy ktoś robił coś takiego w bascomie (jak wogule da rade to obrobić) albo ewentualnie w asemblerze

    jak nie to chociarz jakieś materiały ale przystępnie napisane
  • #18 697778
    h-doc
    Poziom 27  
    W bascomie na pewno da radę, pytani tylko, czy da radę w czasie rzeczywistym. Jeśli nie potrzebujesz zbyt dużej rozdzielczości widmowej, to spróbuj na początku zwykłą DFT, która jest dużo prostsza w implementacji niż FFT, ale też znacznie wolniejsza.
    W asmeblerze znajdziesz również przykłady, ale raczej na procki DSP.
  • #19 698152
    morph13
    Poziom 25  
    Przyznam, że śledzę temat DFT i FFT na forum i jak dotąd nie pojawił się żaden konkretny link lub wstawka programowa dotycząca pełnej implementacji DFT lub FFT w asm dla '51 (dominującego proca na forum), jest to temat "chodliwy" i jeśli ktoś ma jakies osiagnięcia w tej dziedzinie to je skrzrtnie skrywa, poniekąd rozumien to bo takie są prawa rynku :-(.

    Dlatego zbieram się sam, żeby właśnie na bazie plików asm z DSP (320F241) stworzyć coś samemu, no myślę, że jeszcze pomoże mi SIEGER :-)

    Jeszcze raz proponuje przerzucić temat do "DSP i Transmisja" może bardziej się rozwinie.
    Już przeniesiono
    j.
  • #20 703845
    hunterhouse
    Poziom 26  
    nikt chyba nie ma fft w bascomie więc spróbuje sam ale będziecie mi musieli trochę pomóc

    załózmy że mam tablice zmiennych i zczytuje do niej wartości z przetwornika
    ale co ile czasu i ile ich będzie potrzebnych do analizatora widma (7 pasm)

    następnie jakiś prosty wzór do obliczania by się przydał
  • #21 704191
    h-doc
    Poziom 27  
    7 pasm będzie wymagało przeprowadzenia 8-punktowej FFT. Przy czym wtedy masz małą rozdzielczość dla niskich częstotliwości, bo rozdzielczość takiej analizy to fs/n, gdzie fs to częstostliość próbkowania a n liczba punktów fft. Jeśli chodzi o czas potrzebny na przetworzenie danych, to można go obliczyć w przyblizeniu zakładając, że najbardziej czasochłonna jest operacja mnożenia. Dla FFT liczba mnożeń wynosi n/2*log2(n).
  • #22 704336
    hunterhouse
    Poziom 26  
    ale czy mogli byście mi podać konkretne wzory na oblicznie albo jakiś algorytm ale przystępnie napisane
  • #23 705836
    hunterhouse
    Poziom 26  
    no dalej niech mi ktoś pomoże bo ze znalezionych materiałów dużo nie rozumiem

    ????????????????????????????????????????????????????????????????????
  • #24 706286
    h-doc
    Poziom 27  
    słuchaj, jeśli FFT sprawia Ci problem, zacznij od DFT. Tylko jeden prosty wzór:
    gdzie:
    X(k) - k-ta próbka transformaty
    x(n) - n-ta próbka sygnału
    N - całkowita liczba punktów transformaty
  • #25 707485
    hunterhouse
    Poziom 26  
    ten wzór to już coś konkretnego

    oprucz tego znalazłem na stronce http://vitecvitec.webpark.pl/index.htm
    coś takiego nie jest to program ale dobrze obrazuje algorytm DTF

    tylko czy jest poprawny ???

    5 N = 512 // zakładamy długość sygnału wejściowego
    10 DIM SIG [N] // deklaracja tablicy na 512 próbek, N = 512
    20 DIM REX [N/2+1] // deklatacja tablicy ReX o długości N/2+1 = 257
    30 DIM IMX [N/2+1] // deklatacja tablicy ImX o długości N/2+1 = 257
    40 SIG := sygnał wejściowy do analizy // symbolicznie opisane ładowanie tablicy SIG
    // zerowanie tablic ReX i ImX
    50 FOR I:=0 TO (N/2+1)
    60 ReX[I] := 0
    70 ImX[I] := 0
    80 LOOP I
    // wyliczenie właściwej DFT
    90 FOR K:=0 TO (N/2) // dla każdego elementu tablic ReX i ImX
    100 FOR I:=0 TO (N-1) // dla każdej próbki sygnału wejściowego
    110 REX[K[ := REX[K] + SIG[I] * COS(2*pi*I*K / N)
    120 IMX[K[ := IMX[K] - SIG[I] * SIN(2*pi*I*K / N) // pamiętamy o minusie
    130 LOOP I // powtórz dla każdej próbki sygnału wejściowego
    140 LOOP K // powtórz dla każdego elementu tablic ReX i ImX

    jeśli jest to dobrze a jak nie to proszę o jego poprawe

    i tak przy okazji to jak z tech dwóch współczyników REX i IMX wyliczyć moc sygnału czżby wystarczyło je dodać
  • #26 707961
    h-doc
    Poziom 27  
    Program wygląda w porządku (pamiętaj, że oblicza on tylko połowę widma i będzie on działać tylko dla sygnałów rzeczytwistych). Ten link, który podałeś jest bardzo interesujący, myślę, że powinieneś postudiować zawarte tam materiały.
    Odnośnie pytań o widmo:
    - amplituda prążka to moduł z liczby zespolonej, która gi reprezentuje
    - amplituda prążka zależy od całkowitej liczby punktów DFT, więc żeby wyniki były zgodne dla różnych wartości N wprowadza się pojęcie znormalizowanej amplitudy (wartość próbki DFT dzielona jest przez N)
    - odnośnie mocy prążka to nie jestem pewien, ale jest to chyba kwadrat modułu prążka (tutaj znowu możemy wprowadzić moc znormalizowaną)

    Teraz wypadałoby się zainteresować takimi rzeczami jak dobór okna czasowego i nakładki. Powodzenia!
  • #27 709456
    hunterhouse
    Poziom 26  
    ok to dobrze że ten algorytm działa

    mam jeszcze inne pytanie czy mogę podzielić widmo na dwie części i każde z nich osobno prubkować i obliczać
    to po to by mieć dużą rozdzielczość dla basów i duże pasmo dla sopranów

    np 2kHz dla Basów 32 próbki to da rozzielczość 60Hz-1kHz co 60Hz
    i 30kHz dla sopraów 32 próbki to da rozdzielczość 1kHz-15kHz co 1 kHz

    i czy DFT dla 32 próbek będzie dość dokładny do analizatora widma
  • #28 709524
    h-doc
    Poziom 27  
    W zasadzie można obliczyć DFT tylko dla części widma - wystarczy, że obliczysz X(k) tylko dla interesujących Cię argumentów k. Można zatem przeprowadzić częściowe DFT dla poszczedólnych fragmentów pasma z różną rozdzielczością.

    rozdzielczość widmowa DFT wynosi fs/N, gdzie fs to częstotliwość próbkkowania.
  • #30 723947
    shg
    Poziom 35  
    A ja coś wiem! I powiem, nie będę szuja :D .
    Transformate Fouriera ma pewne wady - potrzebuje odrobinę pamięci (DFT), albo nawet więcej (FFT, ale zależy od implementacji), poza tym sama procedura do najkrótszych nie należy. W wypadku użycia FFT trzeba obliczać wszystkie składowe widma. Nie trzeba tego robić dla DFT. Ale co gdy sygnał nie "trafi" na przetwarzane pasmo? Nie zobaczysz go na analizatorze.
    Istnieje alternatywa, myślę że takie rozwiązanie powinno zdać egzamin. Nazywa się filtr IIR. Nie wymaga dużo pamięci (3 bajty na każdy kanał + 3 bajty na dane wejściowe), nie trzeba wykonywać wielu obliczeń i co równie ważne - przetwarza sygnał "w locie", tzn. nie musi gromadzić odpowiedniej (np.256) liczby próbek, wystarczą jej tylko 3 próbki. Posiada szczególnie ważną cechę - charakterystyki tych filtrów są identyczne jak filrów analogowych, więc analizator widma będzie nie do odróżnienia od analogowca :D.
    Dorzucam plik PDF z implementacją 10 pasmowego eqalizera stereo na DSP56001 Motoroli. Podał bym tylko link, ale niestety nie wiem skąd go ściągałem, bo to było ze 2 lata temu :( . Dość dokładnie opisana teoria, aczkolwiek nieco skomplikowana. Plik żródłowy w asm, ale jest też programik w BASICU :!: (pojedynczy filtr, nie cały equalizer i takie właśnie filtry należy wykorzystać) Wyliczone tabele współczynników dla filtrów przy częstotliwości próbkowania 44.1kHz (pasma częstotliwości typowe dla equalizerów, czyli także dla analizatorów widma).
    I jeszcze taka rzecz - po przepuszczeniu sygnału przez te filtry należy przetworzyć otrzymany sygnał tak, aby uzyskać jego amplitudę. A jak przemnoży się te sygnały przez odpowiednie wagi, zsumuje i poda na wejście przetwornika C/A to będzie dopiero equalizer :D .
    Załączniki:
REKLAMA