ledzenie częstotliwoci podstawowej metodš McAulay'a-Quatieri'ego

Marek Szczerba, maj 2000

Inicjalizacja

Wczytanie bibliotek

Biblioteka obsługi dwięku

<< Miscellaneous`Audio`<< Statistics`MultiDescriptiveStatistics`

Pakiet przetwarzania sygnału

Needs["SignalProcessing`"]

Akwizycja sygnału

\(SetDirectory["D:\MIDI\"]; 

    snd = ReadSoundfile["mforts.wav", PrintHeader -> True]; \)

Parametry pracy

Parametry definiowane przez użytkownika

Częstotliwoć próbkowania

CzP=22050;

Rozmiar okna czasowego

Okno=512;

Wielkoć zakładki

Zakladka=(1/2);

Współczynnik zeropaddingu

WspZeroP=1;

Zakres częstotliwoci analizy

CzMin=0;
CzMax = (CzP/2) - (CzP/10);

Ograniczenie liczby składowych harmonicznych (∞ - bez ograniczenia)

MaksHarm=10;

Próg szumowy (0 - wyłšczony)

ProgSz=0;

Próg ciszy

PrCs=-∞;

Kryterium cišgłoci przy odchyleniu częstotliwoci w procentach (MQ - 10%)

OdchCz=(5/100);

Zakres częstotliwoci poszukiwania (Subscript[f, 0])

F0Min=100;
F0Maks=2000;

Parametry wyznaczane automatycznie

Połowa okna czasowego

PolOkna = (Okno/2);

PolOknaAn = PolOkna + PolOkna*WspZeroP;

RelDoRek = (PolOkna/PolOknaAn);

Skok między ramkami

Skok = Okno - Okno*Zakladka;

Liczba ramek analizy

Ramki = Floor[(Length[snd] - Okno*Zakladka)/Skok];

Konwersja zakresu analizy na próbki

\(PrMin = Floor[CzMin\/CzP*Okno]; 
PrMax = Ceiling[CzMax\/CzP*Okno]; 
    \)

Konwersja zakresu poszukiwania (Subscript[f, 0]) na próbki

\(F0PrMin = Floor[\(F0Min\/CzP\) Okno]; 

    F0PrMax = Ceiling[\(F0Maks\/CzP\) Okno]; \)

Detekcja maksymalnej wartoci bezwzględnej próbki w sygnale

MaksPróbka = Max[Abs[snd]];

Definicje funkcji przetwarzania sygnału

Okno czasowe

FOkna=Function[
		z,window=N[Table[DiscreteWindow[Hamming, Okno,nok],{nok, 0, Okno - 1}]];
wy = Table[0, {Length[window]}];
		For[xok = 1, xok <= Length[window], xok++, 
                wy[[xok]] = 
                  z[[xok]]*
                    window[[xok]]];
			wy];

Parametry widmowe (widmo zespolone, amplituda, faza)

WidmoZ = Function[wx, Fourier[wx]];

Amplituda = Function[ax, Abs[ax]];

Faza = Function[px, Phase[px]];

\(\(ZeroP = 
      Function[{zx, \ zy}, 
		ToAdd = Table[0, \ {zz, \ Length[zy]}]; 

        		Flatten[Append[zy, \ ToAdd]]]; \)\)

Obliczenie parametrów widmowych

Obliczenie widma

Obliczenie macierzy widma, amplitudy i fazy

Przed rozpoczęciem obliczeń tworzone sš macierze zerowe o wymiarach Ramki x PolOkna. Nie wykorzystano funkcji Append z uwagi na znaczny czas przeprowadzenia operacji na macierzach.
Następnie kolejne ramki sygnału odczytywane sš ze zmiennej snd i mnożone przez okno czasowe (zdefiniowane w funkcji FOkna). Następnie obliczane jest widmo dla każdej ramki, które zapisywane jest w zmiennej TabWidma. TabAmpl i TabFazy zawierajš odpowiednio amplitudy i fazy widma.
Iteracja ze zmiennš j jest wykorzystywana do prawidłowego przepisania widma ramek bieżšcych w macierzy.
Na czas tworzenia programu wprowadzono funkcję Timing wywietlajšcš czas prowadzenia obliczeń. Do usunięcia po ostatecznym opracowaniu.

\(Timing[TabWidm = Table[0, {i, Ramki}, {j, PolOknaAn}]; 
\ 
      For[i = 0, i < Ramki, \(i++\), 
	
        AktRamka = 
          Table[snd[[i*Skok + j]], {j, 
              Okno}]; 
		
        If[20\ Log[10, Max[Abs[AktRamka]]\/MaksPróbka] >= PrCs, 

          			OknRamka = FOkna[AktRamka]; 
			
          ZeroPOkn = ZeroP[WspZeroP, OknRamka]; 
			
          WidmoRamki = WidmoZ[ZeroPOkn]; 
			
          For[j = 1, j <= PolOknaAn, \(j++\), 
				
            TabWidm[[\ i + 1, j]] = 
              WidmoRamki[[j]]]]]]\)

\(TabAmpl = Abs[TabWidm]; 
TabFazy = Phase[TabWidm]; \)

Obliczenie maksimów widma

Dla każdej ramki znajdowane sš maksima lokalne w zakresie częstotliwoci zdefiniowanym zmiennymi PrMin i PrMax. Maksima interpolowane sš funkcjš kwadratowš (funkcja Fit). Tworzona jest tablica wyjciowa Maksima zawierajšca wartoci maksimów i odpowiadajšce im częstotliwoci (w skali dyskretnej) i fazy dla poszczególnych ramek.

\(Maksima = {}; 

    For[i = 1, i <= Ramki, \(i++\), 
	AktRamka = {}; 
		
      For[j = \((PrMin - 1)\)\/RelDoRek, \ j <= \((PrMax + 1)\)\/RelDoRek, 
          \(j++\), 
			
          If[TabAmpl[[i, j]] > 
                TabAmpl[[i, j - 1]] && 
              TabAmpl[[i, j]] > 
                TabAmpl[[i, j + 1]], 

            				DoInt = 
              Table[TabAmpl[[i, k]], {k, 
                  j - 1, j + 1}]; 
				
            MaksInt = Fit[DoInt, {1, x, x\^2}, x]; 
				
            Maks = FindMinimum[\(-MaksInt\), {x, 0}]; 
				
            AppendTo[
              AktRamka, {\(-Maks[[1]]\), 
                \ j - 3 + x /. 
                  Maks[[2]], \ 
                TabFazy[[i, j]]}]; ]]
	
        AppendTo[Maksima, AktRamka]; 
]; 
Int2Maksima = Maksima; \)

Redukcja macierzy zgodnie z ograniczeniem liczby składowych i progiem szumowym

Macierze zawierajšce dane o maksimach widmowych w poszczególnych ramkach (interpolacja wielomianem) sš analizowane i  usuwane sš wartoci o amplitudzie poniżej progu szumowego Progsz. Następnie, jeli liczba składowych harmonicznych w ramce przekracza wartoć MaksHarm, usuwane sš wartoci o najniższej amplitudzie.

\(Parametry = Int2Maksima; 

    For[i = 1, i <= Ramki, \(i++\), 
	
      For[j = 1, 
        j <= Length[Parametry[[i]]], , 

        		If[Parametry[[i, j, 1]] < 
            ProgSz, Parametry = \ Delete[Parametry, {i, j}], \(j++\)]]; 
	
      While[Length[Parametry[[i]]] > 
          MaksHarm, 
		MinWart = ∞; 
		
        For[j = 1, 
          j <= Length[Parametry[[i]]], 
          \(j++\), 
					
          If[Parametry[[i, j, 1]] < 
              MinWart, 
						
            MinWart = 
              Parametry[[i, j, 1]]; 

            						MinLok = j; ]]; 
		
        Parametry = Delete[Parametry, {i, MinLok}]; ]
]\)

Następnym krokiem jest utworzenie wektorów zawierajšcych informacje o cišgłoci pršżków w dziedzinie widma. Na podstawie parametrów procentowych cišgłoci odnajdywane sš takie maksima, które spełniajš warunki amplitudowe i częstotliwociowe.

Znalezienie maksymalnej amplitudy widma (do normalizacji)

\(MaksAmp = 0; 

    For[i = 1, i <= Ramki, \(i++\), 
	
      \(For[j = 1, 
        j <= Length[Parametry[[i]]], 
        \(j++\), 
		
        \(If[Parametry[[i, j, 1]] > 
            MaksAmp, \ 
          MaksAmp = 
            Parametry[[i, j, 1]]]; \)]; 
      \)]; \)

Utworzenie tabeli linii łšczšcych pršżki

\(Linie = {}; 

    For[i = 1, i < Length[Parametry], \(i++\), 
	AppendTo[Linie, {}]; 
	
      For[j = 1, 
        j <= Length[Parametry[[i]]], 
        \(j++\), 
		
        A = {Parametry[[i, j, 1]], 
            Parametry[[i, j, 2]]}; 
		
        MinOdchAmp = ∞; 
		
        For[k = 1, 
          k <= Length[
              Parametry[[i + 1]]], 
          \(k++\), 
			
          B = {Parametry[[i + 1, k, 1]], 
              Parametry[[i + 1, k, 2]]}; 
          
			If[
            Abs[A[[1]] - 
                  B[[1]]] > MinOdchAmp, 
            Continue[]]; 
			
          If[Abs[A[[2]] - 
                  B[[2]]] > 
              OdchCz \( A[[2]] + 
                    B[[2]]\)\/2, 
            Continue[]]; 
			
          MinOdchAmp = 
            Abs[A[[1]] - 
                B[[1]]]; 
			
          Ampli = \(A[[1]] + 
                B[[1]]\)\/\(2  MaksAmp
                \); 
			
          Punkt = {A[[2]], 
              B[[2]], Ampli, 
              MinOdchAmp}; ]; 
		
        AppendTo[Linie[[i]], Punkt]; ]]; 
    
For[i = 1, i < Length[Parametry], \(i++\), 
	
      For[j = 1, 
        j <= Length[Parametry[[i + 1]]], 
        \(j++\), 
		
        A = {Parametry[[i + 1, j, 1]], 
            Parametry[[i + 1, j, 2]]}; 

        		MinOdchAmp = ∞; 
		
        For[k = 1, 
          k <= Length[Parametry[[i]]], 
          \(k++\), 
			
          B = {Parametry[[i, k, 1]], 
              Parametry[[i, k, 2]]}; 

          			If[
            Abs[A[[1]] - 
                  B[[1]]] > MinOdchAmp, 
            Continue[]]; 
			
          If[Abs[A[[2]] - 
                  B[[2]]] > 
              OdchCz \( A[[2]] + 
                    B[[2]]\)\/2, 
            Continue[]]; 
			
          MinOdchAmp = 
            Abs[A[[1]] - 
                B[[1]]]; 
			
          Ampli = \(A[[1]] + 
                B[[1]]\)\/\(2  MaksAmp
                \); 
			
          Punkt = {B[[2]], 
              A[[2]], Ampli, 
              MinOdchAmp}; ]; 
		
        AppendTo[Linie[[i]], Punkt]; ]]
      \)

Usunięcie "rozgałęzień" linii

\(LinieB = {}; 

    For[i = 1, i <= Length[Linie], \(i++\), 
	
      AktRamka = Linie[[i]]; 
	
      AktLinie = {}; 
	
      For[j = 1, j <= Length[AktRamka], \(j++\), 
		
        MinOdchAmp = ∞; 
		
        For[k = 1, k <= Length[AktRamka], \(k++\), 
			
          \(If[AktRamka[[j, 1]] == 
                AktRamka[[k, 1]] && 
              AktRamka[[k, 4]] < 
                MinOdchAmp, 
				
            MinOdchAmp = 
              AktRamka[[k, 4]]; 

            				AktLinia = 
              AktRamka[[k]]; ]; \)

          			 (*\ 
            If[AktRamka[[k, 1]] > 
                AktRamka[[j, 1]], 
              Break[]]*) \ ]; 
		AppendTo[AktLinie, AktLinia]; ]; 
	
      AktRamka = AktLinie; 
	AktLinie = {}; 
	
      For[j = 1, j <= Length[AktRamka], \(j++\), 
		
        MinOdchAmp = ∞; 
		
        For[k = 1, k <= Length[AktRamka], \(k++\), 
			
          \(If[AktRamka[[j, 2]] == 
                AktRamka[[k, 2]] && 
              AktRamka[[k, 4]] < 
                MinOdchAmp, 
				
            MinOdchAmp = 
              AktRamka[[k, 4]]; 

            				AktLinia = 
              AktRamka[[k]]; ]; \)

          			 (*\ 
            If[AktRamka[[k, 2]] > 
                AktRamka[[j, 2]], 
              Break[]]\ *) \ ]; 
		AppendTo[AktLinie, AktLinia]; ]; 

      AppendTo[LinieB, AktLinie]; ]\)

Uzupełnienie brakujšcych, pojedynczych pršżków

\(Punkty = {}; 

    For[i = 1, i <= Ramki, \(i++\), 
	AppendTo[Punkty, {}]; 
	
      For[j = 1, j = Ramki, \(j++\), 
		Jest = False; 
		
        For[k = 1, 
          k <= Length[LinieB[[i]]], 
          \(k++\), 
			
          \(If[Parametry[[i, j, 1]] == 
              LinieB[[i, k, 2]], 

            				Jest = True; 
				Break[]; ]; \)]; 
		
        If[Jest == False, \ 
          AppendTo[Punkty[[i]], 
            \ {Parametry[[i, j, 2]], 
              Parametry[[i, j, 1]]}]]]]
      \)

Konwersja próbki - Hz

\(LinieHz = LinieB; 
Wsp = CzP\/\(2  PolOknaAn\); 

    For[i = 1, \ i <= Length[LinieHz], \(i++\), 
	
      For[j = 1, 
        j <= Length[LinieHz[[i]]], 
        \(j++\), 
		
        LinieHz[[i, j, 1]] *= Wsp; 
		
        LinieHz[[i, j, 2]] *= Wsp; ]]
      \)\(PunktyHz = Punkty; 

    For[i = 1, i <= Length[Punkty], \(i++\), 
	
      For[j = 1, 
        j <= Length[Punkty[[i]]], 
        \(j++\), 
				
        \(PunktyHz[[i, j, 1]] *= Wsp; 
        \)]]		\)

Przygotowanie tabeli zawierajšcej wykresy linii

\(\(Wykresy = {}; 

    For[i = 1, \ i <= Length[LinieHz], \(i++\), 
	
      For[j = 1, 
        j <= Length[LinieHz[[i]]], 
        \(j++\), 
		
        AppendTo[Wykresy, 
          Graphics[{
              RGBColor[
                1 + 0  LinieHz[[i, j, 3
                      ]], 
                1 - LinieHz[[i, j, 3]], 
                0.5 - 0.5  LinieHz[[i, j, 3
                      ]]], 
              Thickness[
                0.004*LinieHz[[i, j, 3]]^
                    \(1\/5\)], 
              Line[{{i, 
                    LinieHz[[i, j, 1
                      ]]}, {i + 1, 
                    LinieHz[[i, j, 2
                      ]]}}]}]]]]; 

    For[j = 1, 
      j <= Length[PunktyHz[[i]]], 
      \(j++\), 
		
      AppendTo[Wykresy, 
        Graphics[{
            RGBColor[
              1 + 0  PunktyHz[[i, j, 2
                    ]], 
              1 - PunktyHz[[i, j, 2]], 
              0.5 - 0.5  PunktyHz[[i, j, 2
                    ]]], 
            Thickness[
              0.004*PunktyHz[[i, j, 2]]^
                  \(1\/5\)], 
            Line[{{i - 1\/2, 
                  PunktyHz[[i, j, 1]]}, {
                  i + 1\/2, 
                  PunktyHz[[i, j, 1
                    ]]}}]}]]]; 	\)\)

Show[Wykresy,  Frame -> True];

Detekcja częstotliwoci podstawowej

Definicje

Definicja funkcji sinc

\(\(sinc = Function[fx, 
		Sin[Okno  fx\/2]\/\(Okno\ Sin[fx\/2]\)]; 
    \)\)

Definicja funkcji odległoci

\(\(Dist = 
      Function[fx, 
		If[Abs[fx] > \(2  π\)\/Okno, 0, sinc[fx]]]; 
    \)\)

Definicja obwiedni widma sygnału

\(\(ObwW = 
      Function[{fA, fω}, 
		
        fdata = Table[{fω[[fx]], 
              fA[[fx]]}, {fx, 
              Length[fA]}]; 
		
        For[fi = 1, fi <= CzP, fi *= N[\@\(2\^2\)\%12], 
          AppendTo[fdata, {fi, 0}]]; 
		fdata = Sort[fdata]; 
		
        fint = Interpolation[fdata, InterpolationOrder -> 1]; 
		fint

        	]; \)\)

Definicja funkcji optymalizacji

\(\(OptF = 
      Function[{ω0, fA, fω}, 
		
        fR = Floor[CzP\/\(2\ ω0\)]; 
		fAE = ObwW[fA, fω]; 

        		Wy = 
          ∑\+\(fk = 1\)\%fR\((
              \((fAE[fk\ ω0])\)\ 
                \((∑\+\(fl = 1\)\%MaksHarm\((
                      fA[[fl]]\ 
                          Dist[fω[[fl
                                ]] - fk\ ω0] - 
                        fAE[\ ω0]\/2)\))\))\); 
		Wy]; \)\)

Obliczenie funkcji optymalizacji dla poszczególnych wartoci ω0

\(\(OptAll = 
      Function[dane, 
		Wyjcie = {}; 
		
        fA = Table[
            dane[[fx, 1]], {fx, 
              Length[dane]}]; 
		
        fω = 
          Table[dane[[fx, 2]] 
              CzP\/\(2\ PolOknaAn\), {fx, Length[dane]}]; 
		
        For[fω0 = F0Min, fω0 <= F0Maks, 
          fω0 = N[fω0 \@ 2\%12], 
			
          Opt = N[OptF[fω0, fA, fω]]; 
			
          AppendTo[Wyjcie, {fω0, Opt}]]; 
		Wyjcie]; \)\)

Definicja funkcji błędu

\(\(BłšdProc = 
      Function[dane, 
		OK = 0; 
		fj = 1; 
		
        Błędzik = {}; 
		Liczba = Length[dane]; 
		
        For[fi = 1, fi <= Liczba, \(fi++\), 
			
          While[\(\((fi - 1)\) Skok + 1\)\/CzP > 
                Fraza[[fj, 2]] && 
              fj < Length[Fraza], \(fj++\)]; 
			
          GrDolna = 
            Fraza[[fj, 1]] 
              \@\(2\^\(-1\)\)\%12; 
			
          GrGórna = 
            Fraza[[fj, 1]] \@2\%12; 

          			AppendTo[Błędzik, 
            Fraza[[fj, 1]] - 
              dane[[fi]]]; 
			
          If[GrDolna < dane[[fi]] && 
              dane[[fi]] < 
                GrGórna, \(OK++\)]; 
			
          If[dane[[fi]] == 0, 
            \(Liczba--\)]]; 
		100  OK\/Liczba // N]; \)\)

Obliczenia

Obliczenie funkcji optymalizacji dla wszystkich ramek

\(FS = {}; 

    For[i = 1, i <= Ramki, \(i++\), 
	
      If[Length[Parametry[[i]]] > 0, 

        		Opty = 
          OptAll[Parametry[[i]]]; 
		
        MinOpt = {0, 0}; 
		
        For[j = 1, j <= Length[Opty], \(j++\), 
			
          If[Opty[[j, 2]] < 
              MinOpt[[1]], 
            MinOpt = {Opty[[j, 2]], 
                j}]]; 
		
        AppendTo[FS, 
          Opty[[MinOpt[[2
              ]], 1]]]]]; \)

Konwersja czasu

InS = 
      Table[{(((i - 1) Skok + 1)/CzP), 
          FS[[i]]}, {i, Length[FS]}];

a = 
      ListPlot[InS, PlotJoined -> True, PlotRange -> All, Frame -> True, 
        PlotStyle -> RGBColor[0, 0, 1]];

Show[a, Graph];

Obliczenie błędu detekcji

Błšd = BłšdProc[FS]

\(\(ListPlot[Błędzik, PlotJoined -> True, 
      PlotStyle -> RGBColor[0, 0, 1]]; \)\)

Wartoć rednia miary błędu

\(Mean[Abs[Błędzik]]\)

Rekonstrukcja sygnału

Rekonstrukcja sygnału odbywa się na podstawie wygenerowania poszczególnych ramek na podstawie parametrów zapisanych w tablicy Parametry

\(Timing[
	PopOkno = Table[0, {j, Okno}]; 
Wy = {}; 

      For[i = 1, i <= Ramki, \(i++\), 
		
        L = Length[Parametry[[i]]]; 

        		A = Table[
            Parametry[[i, x, 1]], {x, 
              L}]; 
		
        ω = 
          Table[Parametry[[i, x, 2]]*
              RelDoRek, {x, L}]; 
		
        ϕ = 
          Table[Parametry[[i, x, 3]], {
              x, L}]; 
		δ = N[\(2  π\ \)\/Okno]; 
		
        Rek = Table[
            ∑\+\(n = 1\)\%L\((
                A[[n]] 
                  Cos[\((\ 
                      j\ ω[[n]] + 
                        ϕ[[n]])\)])
                \), {j, δ, 2  π, δ}]; 
		
        RekOkno = FOkna[Rek]; 
		
        Komb = Table[
            RekOkno[[j]] + 
              PopOkno[[j]], {j, Skok}]; 
        
		PopOkno = 
          Table[RekOkno[[j]] + 
              PopOkno[[j]], {j, 
              Skok + 1, Okno - Skok}]; 
		
        AppendTo[PopOkno, 
          Flatten[Table[
              RekOkno[[j]], {j, 
                Okno - Skok + 1, Okno}]]]; 
		PopOkno = Flatten[PopOkno]; 
        
	AppendTo[Wy, Flatten[Komb]]; ]; 

      AppendTo[Wy, 
        Table[RekOkno[[j]], {j, 
            Okno - Skok, Okno}]]; 
Wy = Flatten[Wy]; ]\)

ListPlay[Wy, PlayRange -> All, SampleDepth -> 16,  
      SampleRate -> 22050];

Oryginał = 
      Table[snd[[j]], {j, Ramki*Skok}];

ListPlay[Oryginał, PlayRange -> All, SampleDepth -> 16, 
      SampleRate -> 22050];