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

symulacyjny model filtru polifazwoego

cepa 16 Gru 2004 01:52 67035 164
  • #151 16 Gru 2004 01:52
    And!
    Admin grupy Projektowanie

    U mnie jest na dynamicznych.

    Może to zabrzmi śmiesznie ale zadeklaruj na początek tablicę o rozmiarze (2*n)+2
    gdzie n to ilość danych jakie wpiszesz do tablicy. (reszte wyzeruj)

    U mnie notorycznie wychodził poza zakres.

    Zadeklarowałem większą tablicę. :)

  • Mitronik
  • #152 18 Gru 2004 16:35
    And!
    Admin grupy Projektowanie

    Co do samego wyświetlania to używałem takiego sposobu:

    Form1->Canvas->Pixels[x][y] =TColor(kolor) ;

    Jest to dość powolny sposób.

    Natomiast po zastosowaniu Timage zaproponowanej przez myrcioch tutaj
    https://www.elektroda.pl/rtvforum/topic206896.html
    Szybkość wyświetlania polepszyła się znacząco.
    Załączam zmodyfikowany program.
    Stary program jest tu:
    https://www.elektroda.pl/rtvforum/topic14925-0-asc-120.html

    duży plik dla nie mających Buildera
    mały dla mających odpowiednie biblioteki

  • #153 23 Sty 2005 21:38
    RaInMan
    Poziom 11  

    czy ktos mi moze powiedziec dlaczego po FFT jak zeruje np. jedna czestotliwosc, nastepnie robie IFFT i zamiast czegos podobnego co mialem na wejsciu otrzymuje na wyjsciu sam szum :\

  • #154 23 Sty 2005 22:44
    shg
    Specjalista techniki cyfrowej

    a jak nie wyzerujesz to jest OK?
    Może być spowodowane tym, że w ramach Powszechnej Akcji Zwiększania Szybkości procedura liczące FFT nie skaluje danych wyjściowych, natomiast IFFT potrzebuje danych "wyskalowanych", stąd mamy niezłe przesterowanie. Przykład:
    Dla N elementowego wektora wejściowego (w) po wykonaniu funkcji
    v:=fft(w)
    trzeba w MathCADzie pomnożyć wektor v przez 1/√N , żeby uzyskać prawidłowe wyniki, w innych implementacjach może być inaczej, np. wynik z transformaty wykonanej za pomocą funkcji z pakietu kiss_fft trzeba mnożyć (o ile dobrze pamiętam) przez 1/N i dopiero takie dane "przepuszczać" przez ifft.

    W mathCADzie są z tej przyczyny zdefiniowane dwie pary transformat: fft()-ifft() i FFT()-IFFT(), pierwsza para działa bez zkalowania, druga ze skalowaniem

    Najlepiej sprawdź, jakie wartości "wypluwa" fft, wtedy wszystko będzie jasne.

    Ciekawostka: wszystkie fft, z jakimi się spotkałem, w wersji operującej na liczbach rzeczywistych wprowadzają niewielki błąd do składowej urojonej (poniżej 1ppm), natomiast wersje operujące na liczbach zespolonych podają dokładne wartości :D

  • Mitronik
  • #155 23 Sty 2005 23:04
    RaInMan
    Poziom 11  

    wydaje mi sie ze to o skalowanie nie chodzi bo jak nic nie zmieniam i po FFT robie od razu IFFT to jest wszystko dobrze, ale wystarczy ze cos po tym FFT zmienie to od razu szum :?

  • #156 24 Sty 2005 00:18
    shg
    Specjalista techniki cyfrowej

    Hmm, czyli coś dziwnego. możesz podrzucić jakiś przykładowy sygnał?
    Jaki masz rozmiar transformaty? którą składową wycinasz? A najlepiej źródło, ale samą obróbkę, zakładamy, że fft działa :D .

    Ja zrobiłem taki test na szybko:
    załadowałem plik .wav, przepuściłem cały plik naraz przez fft (N=65536) i wyciąłem składniki widma od 100 do 2000 (czyli tak na oko od 50 do 1000 Hz) i zapisałem spowrotem do wava, plik MathCADa i próbki w załączniku, a jak nie masz MathCADa to na obrazku jest zrzut ekranu.

    właściwie to trochę oszukiwałem (albo i nie), bo nie podzieliłem danych na mniejsze kawałki, tylko zrobiłem transformatę na całym wektorze naraz, ale jak widzać (słychać :D ) szumu brak.

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

    Sprawdziłem jeszcze jedną rzecz:
    zrobiłem sobie sinusa 5.5 okresu na transformatę (N=512) i wyciąłem widmo od połowy w górę, razultat - bez zmian. Zrobiłem to samo z cosinusem i na krańcach przebiegu (mniej więcej po jednym okresie z początku i końca) pojawiły się dość spore zniekształcenia, ale, to też nie może być źródłem szumu, tylko trzasków na granicy przetwarzanych buforów.

  • #157 24 Sty 2005 23:51
    RaInMan
    Poziom 11  

    korzystalem z FFT i IFFT z tego linka http://www.library.cornell.edu/nr/cbookcpdf.html i cos mi sie wydaje ze mimo iz jak zrobimy fft pozniej ifft to jest to samo na wyjsciu co na wejsciu to cos jest z tym fft nie tak bo np. zmieniajac tylko jedna skladowa czestotliwosc z np. fft 65536 otrzymujemy na wyjsciu jakies dziwne rzeczy. Sprobuje zmienic ta funkcje fft na jakas inna i zobaczymy.

  • #159 25 Kwi 2005 15:02
    ogr
    Poziom 14  

    Przeczytalem "temat" i mowiac szczerze troche sie pogubilem,
    nie zajmowalem sie dotad przetwarzamiem sygnalow analogowych i
    jestem w tej dziedinie calkiem "zielony".Chcialbym zrobic dekoder
    DTMF (atmega8 ,8Mzh) i nie wiem jak sie do tego zabrac:liczyc DFT ,FFT,
    a moze jest inny sposob,prosze o jaka rade (moze przyklad w C)
    Pozdrawiam.

  • #160 25 Kwi 2005 18:36
    And!
    Admin grupy Projektowanie
  • #161 25 Kwi 2005 19:52
    shg
    Specjalista techniki cyfrowej

    http://www.dattalo.com/technical/theory/dtmf.html

    DFT będzie szybsza od FFT :) a to dla tego, że nie potrzebujemy całego widma, a tylko prążki ośmiu częstotliwości charakterystycznych dla DTMF.

    Jakby to na '51 zrobić to wychodzi taniej niż specjalizowany dekoder :)

    Na bazie tego można by też zrobić koder DTMF, tylko trzeba by użyć IDFT. Już na forum pisałem o tym, tylko w jeszcze bardziej uproszczonej formie - dwa osobne generatory (programowe) z wyjściami sumowanymi na rezystorach.

  • #162 26 Kwi 2005 08:07
    ogr
    Poziom 14  

    Specjalizowany dekoder nie wchodzi w gre,nie ze wzgledu na koszty,
    poprostu chce to zrobic programowo i miec z tego duzo satysfakcji.
    Moglbym kupic np MT8870 lub caly KIT ,ale nie o to chodzi.
    Troche poszukalem w sieci i mysle ze DFT wystarczy aby to pewnie dziwlalo.Czy probkowanie z czestotliwoscia 4kHz i liczenie DFT z 8 probek
    wystarczy?

  • #163 08 Cze 2005 21:10
    dykan
    Poziom 9  

    witam was !!! czy zna sie ktoś z Was na programie Gnu Octave?? potrzebuje w tym programie napsiac krótki programik obrazujący istotę alisingu oraz likwidację go poprzez filtr dolnoprzepustowy czyli antyaliasing;)) porsze was o pomoc Pomózcie mi bardzo !! dziekuję;))

  • #164 08 Sty 2006 20:08
    stopi
    Poziom 28  

    [quote="And!"][code]void FFT(float dane[], unsigned int dlugosc){
    unsigned int n,mmax,m,j,istep,i;
    float wtemp,wr,wpr,wpi,wi,theta;
    float tempr,tempi,temp;

    n=dlugosc << 1;
    for (i=1,j=1;i<n;i+=2) {
    if (j > i) {
    //SWAP
    temp=dane[j];
    dane[j]=dane[i];
    dane[i]=temp;
    //SWAP
    temp=dane[j+1];
    dane[j+1]=dane[i+1];
    dane[i+1]=temp;
    }

    m=dlugosc;
    while (m >= 2 && j > m) {
    j -= m;
    m >>= 1;
    }
    j += m;
    }
    mmax=2;
    while (n > mmax) {
    istep=mmax << 1;
    theta=(6.28318530717959/mmax); //2PI
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for (m=1;m<mmax;m+=2) {
    for (i=m;i<=n;i+=istep) {
    j=i+mmax;
    tempr=wr*dane[j]-wi*dane[j+1];
    tempi=wr*dane[j+1]+wi*dane[j];
    dane[j]=dane[i]-tempr;
    dane[j+1]=dane[i+1]-tempi;
    dane[i] += tempr;
    dane[i+1] += tempi;
    }
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
    }

    //liczenie modułów
    n=dlugosc/2;
    for (j=1,i=1;j<n;j++,i+=2) {
    temp=sqrt(abs((dane[i]*dane[i])+(dane[i+1]*dane[i+1])));
    dane[j]=temp;
    }

    }
    [/code]

    Może ktoś przetestuje ten kod.............


    duży plik dla nie mających Buildera
    mały dla mających odpowiednie biblioteki[/quote]

    Słuchajcie,ten kod jest ok,testowałem go,tylko...on wymaga tablicy iczb zespolonych wejściowych indeksowanych od 1!
    Tak więc pierwsza dana wejściowa tablicy to nie przykładowo we[0] tylko we[1].Dlatego aby użyć powyższej funkcji musimy przesunąć nasz bufor wejściowy o 1 "do przodu" .Przykładowo dla N=256 funkcja przlicza próbki nie od 0..511 tylko od 1..512.
    Napisałem własną FFT która liczy w 100% dobrze i indeksuje próbki od 0,czyli tak jak w zasadzie powinno się prawidłowo programować w języku C++,bo przecież dla komputera pierwszym indeksem jest zawsze 0.

    Funkcja którą tu podaje wygląda następująco:

    double PI=3.14159265358979323846;
    double test_buf[16];
    double buf[8];

    void fft(double data[],unsigned int N)
    {
    unsigned int n,m,j,step,i,GS,k,a,W_N,x;
    double tmp;
    double WI,WR,theta;
    double temp1,temp2,temp3,temp4,temp1_2,temp3_4,Ur,Ui;

    n=N<<1;




    j=1;
    for(i=1;i<n;i+=2)
    {
    if (j > i)
    {
    //zamien czesc rzeczywista
    tmp=data[j-1];
    data[j-1]=data[i-1];
    data[i-1]=tmp;
    //zamien czesc urojona
    tmp=data[j];
    data[j]=data[i];
    data[i]=tmp;
    }
    m=N;
    while (m >= 2 && j > m)
    {
    j -= m;
    m >>= 1;
    }
    j += m;
    }

    n=N<<1; //n=16 dla 8 liczb zespolonych bo kazda liczba sklada sie z dwoch czesci - real oraz imag

    WR=1.0; //cos(0)=1.0;
    WI=0.0; //sin(0)=0.0;
    theta=2*PI;
    W_N=2;
    step=2;

    x=0;
    a=1;

    GS=n/4; //Group Step = 4
    for(k=(n/4);k>0;k=k/2) //k=4,2,1 (dla N=8 czyli 8 liczb ZESPOLONYCH czyli..16 liczb rzeczywistych)
    {
    for(m=0;m<n;m+=GS)
    {

    for(i=m;i<m+((GS/2)-1);i+=2) //i - I skladnik motylka [Upper]
    { //j - II skladnik motylka [Lower]

    WR=cos((theta*x)/W_N);
    WI=-sin((theta*x)/W_N);

    j=i+step; //odleglosc miedzy skrzydelkami motylka
    temp1 = data[j] * WR;
    temp2 = data[j+1] * WI;
    temp3 = data[j] * WI;
    temp4 = data[j+1] * WR;

    temp1_2 = temp1 - temp2;
    temp3_4 = temp3 + temp4;

    Ur=data[i];
    Ui=data[i+1];

    data[i] = temp1_2 + data[i];
    data[i+1] = temp3_4 + data[i+1];
    data[j] = Ur - temp1_2;
    data[j+1] = Ui - temp3_4;

    x++;
    if( x==a )x=0;

    }
    x=0;
    }
    a=a<<1;
    step=step<<1;
    W_N=W_N<<1;
    GS=GS<<1;
    x=0;
    } //for k....


    //liczenie amplitud
    /* n=16;
    for (i=0;i<n;i+=2)
    {
    tmp=sqrt(abs((data[i]*data[i])+(data[i+1]*data[i+1])));
    data[i]=tmp;
    data[i+1]=0;
    } */
    }

    Oczywiscie bufor wejsciowy dla N=4 bedzie zawierał 4 liczby zespolone,z czego kazda sklada sie z 2ch częsci - urojonej i rzeczywistej,czyli bedzie zawieral 8 indeksów (od 0 do 7).Jeśli bufor wejsciowy jest rzeczywisty (a tak przecież jest praktycznie zawsze,bo probkujemy sygnal rzeczywisty)to po prostu uzupełniamy każdą próbkę częscią urojoną równą 0 i formatujemy bufor tak,aby np.buf[i]="jakas wartosć" a buf[i+1]=0,czyli....
    dla 'i' parzystych gdzie i=2k (gdzie k=0,1,2,3...) buf [i] = "jakas wartosc" ,a dla i=2k+1 (nieparzyste) buf[i]=0.Tak więc faktyczna dana wejsciowa to w postaci zespolonej buf[i] + j*buf[i+1] ,a że buf[i+1] jest rowne 0 to jest ona dana rzeczywistą jak gdyby w zespolonej postaci.W obliczaniach mojej funkcji traktujemy ją oczywisćie jak daną zespoloną czyli bierzemy do obliczen tez część urojoną równą 0.

    Procedurę testowałem na buforze złożonym z 8 zadanych przeze mnie liczb zespolonych.Pisałem to niedawno na procek DSP Texas'a TMS320C6713.To taka "skromna" maszynka wykonująca 1024=punktową FFT w czasie rzeczywistym.Testowałem wyniki pod debugerem,a sprawdzałem z wynikami obliczonymi ręcznie dla tych danych w bardzo dobrej książce pana Richarda Lyons'a o czyfrowym przetwarzaniu sygnałów.Po dwóch tygodniach w końcu wyniki zaczęły przypominać te które miały wyjść :)

    A propo's algorytmu sortowania,to do tego wszystkiego ,przed przeprowadzeniem FFT trzeba dokonywać specyficznego sortowania probek wejsciowych co nazywa się odwróceniem bitowym.
    Dla zainteresowanych prześlę materiały w power point'cie o FFT (radix-2-FFT algorithm) i jej strukturze motylkowej i o tym jak napisać ten kod,czyli jak przeliczać wszystkie motylki itd.Nie jest on faktycznie taki rtudny jak się wydaje.
    Mogę przesłać również cały program z małym przykładowym buforem wejściowym (8-punktowa FFT) i podać jakie powinny być wyniki i możecie sobie sprawdzić,że liczy dobrze.Mój kod nie jest jeszcze zoptymalizowany,bo mam wrażenie,że można go jeszcze uprościć sporo,ale muszę nad tym posiedzieć

  • #165 03 Maj 2007 12:41
    rafmos
    Poziom 13  

    Próbowałem ten algorytm napisany przez stopi niestety nie daje takich samych wyników jak ten wyżej.