Analizuj dźwięk za pomocą szybkiej transformacji Fouriera

Próbuję stworzyć graficzny analizator widma w Pythonie.

Obecnie czytam 1024 bajty 16-bitowego dwukanałowego strumienia audio z częstotliwością próbkowania 44100 Hz i uśredniam amplitudę 2 kanałów razem. Więc teraz mam zestaw 256 podpisanych szortów. Chcę teraz wykonać fft na tej tablicy, używając modułu takiego jak numpy, i użyć wyniku do stworzenia graficznego analizatora widma, który na początku będzie miał tylko 32 słupki.

Przeczytałem artykuły na Wikipedii o szybkiej transformacji Fouriera i dyskretnej transformacji Fouriera, ale nadal nie jestem pewien, co reprezentuje otrzymana tablica. Tak wygląda tablica po wykonaniu operacji fft na mojej tablicy przy użyciu numpy:

   [ -3.37260500e+05 +0.00000000e+00j   7.11787022e+05 +1.70667403e+04j
   4.10040193e+05 +3.28653370e+05j   9.90933073e+04 +1.60555003e+05j
   2.28787050e+05 +3.24141951e+05j   2.09781047e+04 +2.31063376e+05j
  -2.15941453e+05 +1.63773851e+05j  -7.07833051e+04 +1.52467334e+05j
  -1.37440802e+05 +6.28107674e+04j  -7.07536614e+03 +5.55634993e+03j
  -4.31009964e+04 -1.74891657e+05j   1.39384348e+05 +1.95956947e+04j
   1.73613033e+05 +1.16883207e+05j   1.15610357e+05 -2.62619884e+04j
  -2.05469722e+05 +1.71343186e+05j  -1.56779748e+04 +1.51258101e+05j
  -2.08639913e+05 +6.07372799e+04j  -2.90623668e+05 -2.79550838e+05j
  -1.68112214e+05 +4.47877871e+04j  -1.21289916e+03 +1.18397979e+05j
  -1.55779104e+05 +5.06852464e+04j   1.95309737e+05 +1.93876325e+04j
  -2.80400414e+05 +6.90079265e+04j   1.25892113e+04 -1.39293422e+05j
   3.10709174e+04 -1.35248953e+05j   1.31003438e+05 +1.90799303e+05j...

Zastanawiam się, co dokładnie reprezentują te liczby i jak przeliczyłbym te liczby na procent wysokości dla każdego z 32 słupków. Czy powinienem również uśredniać razem 2 kanały?

Tablica, którą pokazujesz, to współczynniki transformacji Fouriera sygnału audio. Współczynniki te można wykorzystać do uzyskania zawartości częstotliwości dźwięku. FFT jest zdefiniowany dla funkcji wejściowych o wartościach zespolonych, więc otrzymane współczynniki będą liczbami urojonymi, nawet jeśli wszystkie dane wejściowe są wartościami rzeczywistymi. Aby uzyskać ilość mocy na każdej częstotliwości, musisz obliczyć wielkość współczynnika FFT dla każdej częstotliwości. To nie jest tylko rzeczywisty składnik współczynnika, musisz obliczyć pierwiastek kwadratowy z sumy kwadratu jego rzeczywistych i urojonych składników. To znaczy, jeśli twój współczynnik wynosi a + b * j, to jego wielkość wynosi sqrt (a ^ 2 + b ^ 2).

Po obliczeniu wielkości każdego współczynnika FFT, musisz dowiedzieć się, do której częstotliwości audio należy każdy współczynnik FFT. Punkt N FFT poda zawartość częstotliwościową twojego sygnału przy N równo rozmieszczonych częstotliwościach, zaczynając od 0. Ponieważ twoja częstotliwość próbkowania wynosi 44100 próbek / sek. a liczba punktów w twojej FFT wynosi 256, twój odstęp częstotliwości wynosi 44100/256 = 172 Hz (w przybliżeniu)

Pierwszym współczynnikiem w Twojej tablicy będzie współczynnik częstotliwości 0. Jest to w zasadzie średni poziom mocy dla wszystkich częstotliwości. Reszta współczynników będzie liczyć od 0 do wielokrotności 172 Hz, aż do 128. W FFT możesz mierzyć częstotliwości tylko do połowy punktów próbkowania. Przeczytaj te linki na temat częstotliwości Nyquista i twierdzenia o próbkowaniu Nyquista-Shannona, jeśli jesteś żarłokiem i musisz wiedzieć, dlaczego, ale podstawowym rezultatem jest to, że twoje niższe częstotliwości będą replikowane lub aliasowane w przedziałach wyższych częstotliwości. Tak więc częstotliwości będą rozpoczynać się od 0, zwiększać się o 172 Hz dla każdego współczynnika aż do współczynnika N / 2, a następnie zmniejszać się o 172 Hz aż do współczynnika N - 1.

To powinno wystarczyć na początek. Jeśli chciałbyś mieć dużo bardziej przystępne wprowadzenie do FFT niż te podane na Wikipedii, możesz spróbować Understanding Digital Signal Processing: 2nd Ed. . To było dla mnie bardzo pomocne.

Więc to właśnie reprezentują te liczby. Konwersję na procent wysokości można przeprowadzić przez skalowanie wielkości każdej składowej częstotliwości przez sumę wielkości wszystkich składowych. Chociaż dałoby to tylko reprezentację względnego rozkładu częstotliwości, a nie rzeczywistą moc dla każdej częstotliwości. Możesz spróbować skalować o maksymalną wielkość możliwą dla składowej częstotliwości, ale nie jestem pewien, czy to się dobrze wyświetli. Najszybszym sposobem znalezienia działającego współczynnika skalowania byłoby eksperymentowanie z głośnymi i cichymi sygnałami audio w celu znalezienia odpowiedniego ustawienia.

Na koniec powinieneś uśredniać oba kanały razem, jeśli chcesz pokazać zawartość częstotliwościową całego sygnału audio jako całości. Miksujesz dźwięk stereo do dźwięku mono i pokazujesz połączone częstotliwości. Jeśli chcesz mieć dwa oddzielne wyświetlacze dla częstotliwości prawej i lewej, będziesz musiał wykonać transformację Fouriera osobno na każdym kanale.

Chociaż ten wątek ma już lata, okazał się bardzo pomocny. Chciałem tylko przekazać swoje uwagi każdemu, kto to znajdzie i próbuje stworzyć coś podobnego.

Jeśli chodzi o podział na słupki, nie należy tego robić, jak sugeruje antti, dzieląc dane równo na podstawie liczby słupków. Najbardziej przydatne byłoby podzielenie danych na części oktawowe, przy czym każda oktawa jest dwukrotnie większa od częstotliwości poprzedniej. (tj. 100 Hz to jedna oktawa powyżej 50 Hz, czyli jedna oktawa powyżej 25 Hz).

W zależności od tego, ile chcesz taktów, cały zakres dzielisz na zakresy 1 / X oktawy. Bazując na zadanej częstotliwości środkowej A na słupku, otrzymujesz górną i dolną granicę słupka z:

upper limit = A * 2 ^ ( 1 / 2X )
lower limit = A / 2 ^ ( 1 / 2X )

Aby obliczyć następną sąsiednią częstotliwość środkową, użyj podobnego obliczenia:

next lower =  A / 2 ^ ( 1 / X )
next higher = A * 2 ^ ( 1 / X )

Następnie uśredniasz dane, które mieszczą się w tych zakresach, aby uzyskać amplitudę dla każdego słupka.

Na przykład: Chcemy podzielić na zakresy 1/3 oktawy i zaczynamy od częstotliwości środkowej 1 kHz.

Upper limit = 1000 * 2 ^ ( 1 / ( 2 * 3 ) ) = 1122.5
Lower limit = 1000 / 2 ^ ( 1 / ( 2 * 3 ) ) =  890.9

Biorąc pod uwagę 44100 Hz i 1024 próbki (43 Hz między każdym punktem danych) powinniśmy uśrednić wartości od 21 do 26. (890,9 / 43 = 20,72 ~ 21 i 1122,5 / 43 = 26,10 ~ 26)

(Paski 1/3 oktawy dają około 30 taktów między ~ 40 Hz a ~ 20 kHz). Jak już możesz się domyślić, gdy będziemy podążać wyżej, będziemy uśredniać większy zakres liczb. Niskie słupki zazwyczaj zawierają tylko 1 lub niewielką liczbę punktów danych. Podczas gdy wyższe słupki mogą być średnią z setek punktów. Powodem jest to, że 86 Hz to oktawa powyżej 43 Hz ... podczas gdy 10086 Hz brzmi prawie tak samo jak 10043 Hz.

masz próbkę, której długość w czasie wynosi 256/44100 = 0,00580499 sekund. Oznacza to, że twoja rozdzielczość częstotliwości wynosi 1 / 0,00580499 = 172 Hz. 256 wartości, które otrzymujesz z Pythona, odpowiada częstotliwościom, w zasadzie, od 86 Hz do 255 * 172 + 86 Hz = 43946 Hz. Liczby, które otrzymujesz, są liczbami zespolonymi (stąd „j” na końcu każdej drugiej liczby).

EDITED: POPRAWIONE BŁĘDNE INFORMACJE

Musisz zamienić liczby zespolone na amplitudę, obliczając sqrt (i ² + j ² ), gdzie i i j są częściami rzeczywistymi i urojonymi.

Jeśli chcesz mieć 32 słupki, powinieneś, o ile rozumiem, wziąć średnią z czterech kolejnych amplitud, uzyskując 256/4 = 32 słupki, jak chcesz.

FFT zwraca N wartości zespolonych, z których można obliczyć module=sqrt(real_part^2+imaginary_part^2). Aby uzyskać wartość dla każdego pasma, należy zsumować moduły dotyczące wszystkich harmonicznych w paśmie. Poniżej możesz zobaczyć przykład dotyczący analizatora widma 10 barów. Kod C musi zostać opakowany, aby uzyskać moduł pyd python.

float *samples_vett;
float *out_filters_vett;
int Nsamples;
float band_power = 0.0;
float harmonic_amplitude=0.0;
int i, out_index;

out_index=0;


for (i = 0; i < Nsamples / 2 + 1; i++)       
        {
            if (i == 1 || i == 2 || i == 4 || i == 8 || i == 17 || i == 33 || i == 66 || i == 132 || i == 264 || i == 511)
            {
                out_filters_vett[out_index] = band_power; 
                band_power = 0; 
                out_index++;  
            }

            harmonic_amplitude = sqrt(pow(ttfr_out_vett[i].r, 2) + pow(ttfr_out_vett[i].i, 2));
            band_power += harmonic_amplitude;

        }

Zaprojektowałem i wykonałem cały 10-ledowy analizator widma słupkowego firmy Python. Zamiast użyć biblioteki nunmpy (zbyt dużej i bezużytecznej, aby uzyskać tylko FFT), utworzono moduł Pythona pyd (tylko 27KB), aby uzyskać FFT i podzielić całe widmo audio na pasma.

Dodatkowo do odczytu wyjściowego dźwięku został stworzony moduł loopback WASapi portaudio pyd. Możesz zobaczyć projekt (schemat blokowy) na obrazku 10BarsSpectrumAnalyzerWithWASapi.jpg

Właśnie dodałem film instruktażowy na moim kanale YouTube: jak zaprojektować i wykonać bardzo inteligentny Python Spectrum Analyzer 10 Led Bar