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

Dokładność FFT, niedokładna metoda?

puciaken 04 Cze 2010 19:30 2462 9
REKLAMA
  • #1 8155343
    puciaken
    Poziom 11  
    Jestem dość początkującym programistą, a zagadnienia związane z FFT są dla mnie po części abstrakcyjne, ale :) staram się rozpoznać wys. dźwięku z użyciem FFT. Kod wygląda mniej więcej tak:
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include "kiss_fftr.h"
    #include "spline2.c"
    
      FILE *f;
      double maxprazekhz, maxprazekval=0, sumhz=0;
      long ilewczytano=0, maxprazek=0, i, maxinterp_i=0;
      int ileprobekfft=0;
      float *x, *y, *y2, *hann, wynik, maxinterp_pr=0, temp;
      
      const double PI=3.141592653589793238462643383279502884197169399375105820974944;
      const unsigned char ilefft=10;
      #define N 16
      #define NFFT 8192
      #define ZASIEG 100
        
      unsigned int iloscprobek;
      signed short int buf=0;
      kiss_fftr_cfg cfg;
      kiss_fft_cfg cfg2;  
      kiss_fft_cpx *cx_out;
      kiss_fft_scalar *cx_in;
      
      x = malloc(sizeof(float)*N);
      y = malloc(sizeof(float)*N);
      y2 = malloc(sizeof(float)*N);  
      for (i=0; i<N; i++)
          x[i] = i;
       
       hann = malloc(sizeof(float)*NFFT);
       for (i=0; i<NFFT; ++i)
       {
           hann[i] = sin((PI*i)/(NFFT-1));
           hann[i] *= hann[i];
       }
    
       f = fopen("440.wav", "rb");
    /*
    OTWARCIE PLIKU, ominięcie nagłówka
    */
          cfg = kiss_fftr_alloc(NFFT,0,NULL,NULL); 
          cx_in = (kiss_fft_scalar*) malloc(NFFT*sizeof(kiss_fft_scalar));
          cx_out = (kiss_fft_cpx*) malloc((NFFT/2+1)*sizeof(kiss_fft_cpx));
          
          while (ilewczytano = fread(&buf, sizeof(buf), 1, f))
          {
                if (ileprobekfft != NFFT)
                {
                    cx_in[ileprobekfft] = buf * hann[ileprobekfft];
                    ++ileprobekfft;
                }
                else
                {
                    if (ilefft > ktorafft)
                    {
                        kiss_fftr( cfg , cx_in , cx_out );
                        ileprobekfft = 0;
                        maxprazekval = 0;
                        maxprazek = 0;
                        maxprazekhz = 0;
                        maxinterp_pr = 0;
                        for (i=0; i<NFFT/2+1; ++i) //szukanie maksimum
                        {
                            temp = hypot(cx_out[i].r, cx_out[i].i);
                            if (temp > maxprazekval)
                            {
                               maxprazekval = temp;
                               maxprazek = i;
                            }
                        }
                        for (i=0; i<N; i++) //wpisywanie wartości do tablicy interpolacji
                            y[i] = hypot(cx_out[i+maxprazek-(N/2)].r, cx_out[i+maxprazek-(N/2)].i);
    
                        spline(x, y, N, 0, 0, y2);  
                        for (i=-ZASIEG; i<ZASIEG; ++i) //dość tandetne szukanie maksimum,
                        { // ale przy dużej wartości ZASIEG i tak wyniki są podobne
                            splint(x, y, y2, N, (N/2)+((float)i/ZASIEG), &wynik);                   
                            if (wynik>maxinterp_pr)
                            {
                               maxinterp_pr = wynik;
                               maxinterp_i = i;
                            }
                        }
                        maxprazekhz = (double)samplerate / NFFT * ((double)maxprazek + ((double)maxinterp_i/ZASIEG));
                        printf(" %LfHz\n", (float)maxprazekhz);
                    }
                    else
                        break;
                 }
          }  
          fclose(w);
      }
      fclose(f);
      free(x);
      free(y);
      free(y2);

    Jestem świadom ilości optymalizacji jakie można wykonać, ale nie to jest teraz ważne. Jak widać robię okno Hanning (a przynajmniej tak mi się wydaje), po czym obliczam FFT, szukam interpolowanego maksimum w otoczeniu największego prążka i obliczam jego częstotliwość. Wszystko byłoby ok, gdyby nie to, że dla sinusoidy 440Hz dostaję wynik 440.2Hz/440.3Hz. Biorąc ten sam plik, zaznaczając 4096 pierwszych sampli i okienko Hanning w Adobe Audition dokładność jest dużo większa (440.03Hz). Co zatem robię źle?
  • REKLAMA
  • #2 8156770
    Xitami
    Poziom 29  
    Uśredniłbym widma z nakładających się fragmentów sygnału
    Dokładność FFT, niedokładna metoda?przykładowy obrazek
    i dopiero w takim uśrednionym widmie szukał maksimum
    należy uwzględnić dynamikę sygnału, prążki mniejsze od pewnego poziomu to tylko szum kwantyzacji, uśrednianie poprawia sprawę (o ile się nie mylę, to N widm zadziała jak dodanie log_2(N) bitów do przetwornika A/D)
    można by robić widma z fragmentu uzupełnionego zerami, okienkowanie przed uzupełnieniem zerami! Dłuższe fft - wyższa rozdzielczość.

    Jeszcze coś, chirp-fft. Gdzieś zeżarło mi notatki, ale to jak lupa, można zobaczyć N/2 prążków z prawie dowolnego wycinka widma.
  • REKLAMA
  • #3 8157257
    puciaken
    Poziom 11  
    To nakładanie się jest dość interesujące. Uśrednianie miałoby polegać na liczeniu śr. arytmetycznej modułów, czy każdej części osobno (rzecz. i urojonej)?
    Na te prążki poniżej pewnego poziomu nie patrzę, więc chyba nie robią mi różnicy. Co jeszcze do uśredniania, nie żebym się czepiał :) ale dlaczego Audition dla FFT 4096 i tej samej ilości próbek daje lepszy wynik? Robi mniejsze uzupełnione zerami i je uśrednia? A może moja interpolacja jest kiepska. Też nie bardzo rozumiem okienkowania przed uzupełnieniem zerami. Przecież na początku bez wczytania próbek już są zera i uzupełnienie polega chyba tylko na tym, żeby np. zamiast wczytywać 4096 próbek wczytać np. 2048.
  • #4 8158244
    Xitami
    Poziom 29  
    Każdy fragment rozpoczyna się w innej chwili w stosunku do głównej składowej, więc faza musi płynąć, znaczy należy uśredniać moduły, wystarczy sumować nawet kwadraty, bo bezwzględna wartość nie ma chyba znaczenia.

    Zasada taka, bierzemy N próbek (w części pokrywające się z poprzednim i/lub następnym przetwarzanym fragmentem) , odejmujemy średnią, mnożymy przez okno (okno "zeruje" 0 i N-1 próbkę) i wołamy FFT dla np. 2N próbek (należy zadbać o to by ogon był zerowy, bo w tablicy mogą być przecie jakieś śmieci)

    Pomyśl jak wygląda sygnał pomnożony przez okno gdy wartość średnia nie jest zerem.
    Dokładność FFT, niedokładna metoda?
  • #5 8159295
    puciaken
    Poziom 11  
    Xitami napisał:
    Każdy fragment rozpoczyna się w innej chwili w stosunku do głównej składowej, więc faza musi płynąć, znaczy należy uśredniać moduły, wystarczy sumować nawet kwadraty, bo bezwzględna wartość nie ma chyba znaczenia.

    chyba pierwiastek? sqrt(im^2 + re^2) ? Wolałbym dla pewności policzyć tradycyjnie moduł. Mam nadzieję, że nie muszę sprawdzać fazy gł. częstotliwości i na tej podstawie dobierać rozmiar overlapu?
    Cytat:

    Zasada taka, bierzemy N próbek (w części pokrywające się z poprzednim i/lub następnym przetwarzanym fragmentem) , odejmujemy średnią, mnożymy przez okno (okno "zeruje" 0 i N-1 próbkę) i wołamy FFT dla np. 2N próbek (należy zadbać o to by ogon był zerowy, bo w tablicy mogą być przecie jakieś śmieci)

    Pomyśl jak wygląda sygnał pomnożony przez okno gdy wartość średnia nie jest zerem.


    Jak zwykle ciężko mi to przyswoić. Rozumiem biorę np 1024 próbek, robię okno (dbając o skrajne punkty, tylko czy okno ma obejmować 1-4096 czy 1-1024?) i fft o rozmiarze np. 4096, potem biorę kolejne 1024 ale tak, że np. 256 pokrywa się z poprzednimi i cykl się powtarza. Na końcu liczę średnią, szukam maksimum i coś wyjdzie.
    Tylko pominąłem "odejmujemy średnią", czego średnią (fft?) i od czego odejmujemy?

    Przepraszam, ale ciężko mi to pojąć :) Liczę na wyrozumiałość.
  • REKLAMA
  • #6 8159484
    Xitami
    Poziom 29  
    moduł jest z pierwiastkiem, a bez pierwiastka mamy kwadrat, :-) rób "tradycyjnie"

    kichamy na fazę, sumować należy moduły

    0. zerujemy widmo średnie
    powtarzamy
    1.należy policzyć wartość średnią sygnału (1024 p.) i odjąć ją od niego, część to poprzedni kawałek ale orginalny (bez odjętej średniej)
    2.mnożymy przez okno, pr. 0 oraz 1023 staną się zerowe
    liczysz tak: hann[i] = sin((PI*i)/(NFFT-1))^2 coś tu nie gra, sprawdź Hann window
    okno ma obejmować tylko sygnał
    3. wołamy FFT(4*N) po zadbaniu o zera na końcu
    4. dodajemy wynik do widma średniego
  • #7 8159942
    puciaken
    Poziom 11  
    Jak zwykle problemy z przyswojeniem, ale:) :
    Xitami napisał:

    powtarzamy
    1.należy policzyć wartość średnią sygnału (1024 p.) i odjąć ją od niego, część to poprzedni kawałek ale orginalny (bez odjętej średniej)

    Coś nie rozumiem z tą średnią. Mam obliczyć średnią wartość próbki w pobranym sygnale i od wszystkich próbek odjąć tą wartość? I "część to poprzedni kawałek ale orginalny (bez odjętej średniej)" oznacza, że od overlapu mam nie odejmować? Chyba coś źle interpretuję :)
    Cytat:

    2.mnożymy przez okno, pr. 0 oraz 1023 staną się zerowe
    liczysz tak: hann[i] = sin((PI*i)/(NFFT-1))^2 coś tu nie gra, sprawdź Hann window
    okno ma obejmować tylko sygnał

    http://en.wikipedia.org/wiki/Hann_function tu jest "mój" wzór. To chyba to samo?
    Cytat:
    3. wołamy FFT(4*N) po zadbaniu o zera na końcu
    4. dodajemy wynik do widma średniego

    Z tym uśrednianiem widma to ma być tak, że liczę widmo tego co mam aktualnie i widmo średnie = [(iledotejporybyłoliczonychfft - 1)*(widmo średnie) + (aktualnie policzone widmo)] / (iledotejporybyłoliczonychfft), czyli właściwie śr. arytm.
    ?
  • #8 8160074
    Xitami
    Poziom 29  
    Twoje okno jest przekozackie, kiedyś nad nim pomyślę, bo bardzom ciekaw jakie widma, duch i zjawy mogą przez nie przejść :-)
    Śpiącym już bardzo, więc mogłem coś pochrzanić więc ostrożnie
    	const PRZES=256, N=1024, nFFT=4096;
    
    	signed short int buf[N+PRZES]
    	double wSrednie[nFFT/2]
    	
    	for i=0..N-1
    		okno[i]=sin(Pi*i/(N-1))  // twój wzór to coś zupełnie innego
    		okno[i]*=okno[i]
    	for i=N..nFFT-1
    		cx_in[i]=0;  //jeżeli FFT nie miesza w cx_in to wystarczy raz na początku
    	for i=0..nFFT/2
    		wSrednie=0;
    	
    	ilewczytano = fread(&buf, sizeof(buf[0]), N+PRZES, f);
    	if ilewczytano < N
    		dupa zbita
    	for i=0..N-1
    		cx_in[i]=buf[i]
    	do
    		srednia=0
    		for i=0..N-1
    			srednia += cx_in[i]
    		srednia /= N
    		for i=0..N-1
    			cx_in[i] = (cx_in[i]-srednia)*okno[i]
    		FFT(nFFT)
    		for i=0..nFFT/2
    			wSrednie[i] += hypot(cx_out[i].r, cx_out[i].i); //zakładam, że to pierwiastek sumy kwadratów
    		
    		for i=0..N-1
    			cx_in[i]=buf[i+PRZES]
    		ilewczytano = fread(&buf[N], sieof(buf[0]), PRZES, f)
    	while ilewczytano==PRZES

    średnia arytmetyczna to poprostu suma, raz jedyny na końcu podzielona przez liczbę sumowań
  • #9 8160855
    puciaken
    Poziom 11  
    Xitami napisał:
    Twoje okno jest przekozackie, kiedyś nad nim pomyślę, bo bardzom ciekaw jakie widma, duch i zjawy mogą przez nie przejść :-)

    Nie bardzo rozumiem. To co napisałeś, co ja napisałem i co jest na wiki, to jest dokładnie to samo. U mnie zamiast N jest NFFT, bo zbierałem tyle próbek ile liczy NFFT. Poza tym nie widzę żadnej różnicy.
    Cytat:

    średnia arytmetyczna to poprostu suma, raz jedyny na końcu podzielona przez liczbę sumowań

    Wiem, wiem, coś miałem zaćmienie z liczeniem tego na piechotę.

    Jutro postaram się to u siebie zastosować, teraz nie mam czasu. Dzięki za pomoc :)

    EDIT: Nie wiem o co Ci dalej chodzi z funkcją okna, sprawdziłem wzór z cos i wychodzi dokładnie to samo. A co do programu, wyniki wychodzą 441Hz-442Hz w zależności jaki overlap się ustawi (i nie ma tu żadnej proporcjonalności).
  • REKLAMA
  • #10 8198395
    puciaken
    Poziom 11  
    pomoże ktoś jeszcze?
REKLAMA