Detekcja częstotliwoci podstawowej przy pomocy funkcji autokorelacji

Š 1999, 2000, Marek Szczerba

Inicjalizacja

Wczytanie bibliotek

Biblioteka obsługi dwięku

<< Miscellaneous`Audio`

Biblioteki statystyczne

<< Statistics`DescriptiveStatistics`

Akwizycja sygnału

SetDirectory["D:\\Mathematica\\Data\\"];

\(snd = ReadSoundfile["oboe1.wav"]; 
AKFName = "AutoK.dat"; \)

Parametry definiowane przez użytkownika

Odstęp w próbkach pomiędzy poszczególnymi punktami analizy

Skok=100;

Częstotliwoć próbkowania

CzP=22050;

Granice pasma analizy

CzMin=100;
CzMax=2000;

Wartoć parametru tłumienia niskich częstotliwoci

Tłumik=0.8;

Parametry obliczane automatycznie

Konwersja pasma analizy

\(OkrMin = 1;  (*Ceiling[CzP\/CzMax]; *) 
OkrMax = Floor[CzP\/CzMin]; 
    \)

Rzeczywista częstotliwoć podstawowa dwięku

RzCzPodst = N[444 (( (2^Wys))^(1/12))];

Liczba punktów analizy

\(LRamek = 0; 

    For[j = 1, j <= Length[snd] - 4  OkrMax, j += Skok, \(LRamek++\)]; \)

Definicje funkcji przetwarzania sygnału

Definicja funkcji autokorelacji

\(\(AutoCorr = 
	
      Function[dane, 
		fL = Length[dane]\  - \ OkrMax; 
		
        Wynik = Table[0, \ {fi, OkrMax}]; 
		
        init = ∑\+\(fn = 1\)\%fL Abs[
              dane[[fn]]]; 
		
        For[fi = OkrMin, \ fi\  <= \ OkrMax, \ \(fi++\), 
			
          ρ = 
            ∑\+\(fn = 1\)\%fL\((
                dane[[fn]] 
                  dane[[fn\  + \ fi]])
                \); 
			
          Wynik[[fi]] = 
            N[ρ\/init]]; 
\ \ \ \ \ \ \ \ Wynik]; \)\)

Pętla funkcji autokorelacji

\(\(OblAutoK = 
	
      Function[{dzw, lr}, 
		stream = OpenWrite[AKFName]; 
		Numer = 0; 
        
		For[fk = 1, fk <= lr*Skok, fk += Skok, 
			\(Numer++\); 

          			Ramka = 
            Table[dzw[[fj]], {fj, fk, 
                fk + 2\ OkrMax - 1}]; 
			AKRamka\  = \ AutoCorr[Ramka]; 
          
			For[fm = 1, fm < OkrMin, \(fm++\), 
				
            AKRamka[[fm]] = 0]; 
			
          Write[stream, AKRamka]]; 
		Close[stream]; ]; \)\)

Odczytanie wartoci funkcji autokorelacji z pliku

\(\(ReadAK = 
      Function[lr, 
		stream = OpenRead[AKFName]; 
		WyP = {}; 
		
        For[fk = 1, fk <= lr\ Skok, fk += Skok, 
				
          AppendTo[WyP, Read[stream, Expression]]]; 
		Close[stream]; 

        		WyP]; \)\)

Poszukiwanie maksimum funkcji autokorelacji

\(\(MaksAutoK = 
	
      Function[dane, 
		Wynik = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), \ 
			Maks = 0; 
			
          Lok = 0; 
			
          For[fj = 1, 
            fj <= Length[dane[[fi]]], 
            \(fj++\), 
				
            \(If[dane[[fi, fj]] > Maks, 
              
					
              Maks = dane[[fi, fj]]; 

              					Lok = fj; ]; \)]; 
			
          AppendTo[Wynik, {Lok, Maks}]; ]; 
		Wynik]; \)\)

Funkcja tłumienia składowych o niskiej częstotliwoci

\(\(Tłumienie = 
      Function[{dane, k}, 
		Wynik = dane; 
		M = OkrMax - OkrMin; 

        		For[fn = 1, fn <= Length[dane], \(fn++\), 
			
          For[fm = 1, 
            fm <= Length[dane[[fn]]], 
            \(fm++\), 
				
            \(Wynik[[fn, fm]] = 
              dane[[fn, fm]]*
                \((1 - \(k\ fm\)\/M)\); \)]]; 
		Wynik]; \)\)

Wyznaczenie pasm półtonowych

\(\(GetBands = 
      Function[{fin, fmin, fmax}, 
		Wy = {fin}; 
		fx = fin; 
		
        While[fx > fmin, 
			fx *= \@\(2\^\(-1\)\)\%12; 
			
          AppendTo[Wy, N[fx]]]; 
		fx = fin; 
		
        While[fx < fmax, 
			fx *= \@2\%12; 
			AppendTo[Wy, N[fx]]]; 
        
		Sort[Wy]]; \)\)

Wyznaczenie maksymalnych wartoci funkcji autokorelacji w pasmach

\(\(GetMaksBands = 
      Function[{fdata, fbands}, 
		
        Wy = Table[Table[0, {Length[fbands]}], {Length[fdata]}]; 
		
        For[fi = 1, fi <= Length[fdata], \(fi++\), 
			
          For[fj = 1, fj <= Length[fbands], \(fj++\), 
				
            TmpTable = {}; 
				
            fkmax = Ceiling[
                CzP\/\(fbands[[fj]] 
                    \@\(2\^\(-1\)\)\%24\)]; 
				
            If[fkmax > OkrMax, fkmax = OkrMax]; 
				
            For[fk = 
                Floor[CzP
                    \/\(fbands[[fj]] 
                      \@2\%24\)], fk <= fkmax, \(fk++\), 
					
              AppendTo[TmpTable, 
                fdata[[fi, fk]]]]; 

            				Wy[[fi, fj]] = 
              Max[TmpTable]]]; 
		Wy]; \)\)

Pomiar poprawnoci detekcji

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

Obliczenia

Obliczenie funkcji autokorelacji dla poszczególnych ramek sygnału

OblAutoK[snd, LRamek];

Odczytanie wartoci funkcji autokorelacji z pliku

AutoK = ReadAK[LRamek];

Tłumienie wartoci funkcji autokorelacji dla niższych częstotliwoci

TłAutoK = Tłumienie[AutoK, Tłumik]/10000;

\(\(For[i = 1, i <= Length[snd], \(i++\), 
	
      If[snd[[i]] == 0, 
        snd[[i]] = 1]]; \)\)

Segmentacja sygnału

sft = 
      Table[Max[TłAutoK[[j]]]*
          (1 + (1/(10  OkrMax)) 
                (∑_(k = 1)^OkrMax) Log[
                    Abs[snd[[k + 
                          (j - 1) Skok]]]]), {j, 
          Length[TłAutoK]}];

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

\(sfa = {}; 
PrevMaxLoc = {1, 1}; 

    PrevTable = 
      Table[0, {
          Length[TłAutoK[[1
              ]]]}]; 

    For[i = 1, i <= Length[TłAutoK], \(i++\), 
	MaxLoc = {0, 0}; 

      	MaxVal = {
          Max[TłAutoK[[i]]], 0}; 
      
	For[j = 1, 
        j <= Length[
            TłAutoK[[i]]], 
        \(j++\), 
		
        If[TłAutoK[[i, j]] == 
            MaxVal[[1]], 
          MaxLoc[[1]] = j]; 
		
        If[TłAutoK[[i, j]] > 
              MaxVal[[2]] && 
            TłAutoK[[i, j]] < 
              MaxVal[[1]], 
          MaxVal[[2]] = 
            TłAutoK[[i, j]]; 
          MaxLoc[[2]] = j]]; 
	
      sfal = \((
              TłAutoK[[i, 
                MaxLoc[[1
                  ]]]])\)*
            \((PrevTable
                [[MaxLoc[[1
                  ]]]])\) + 
			
          \((TłAutoK[[i, 
                MaxLoc[[2
                  ]]]])\)*
            \((PrevTable
                [[MaxLoc[[2
                  ]]]])\) + 
			
          \((TłAutoK[[i, 
                PrevMaxLoc[[1
                  ]]]])\)*
            \((PrevTable
                [[PrevMaxLoc[[1
                  ]]]])\) + 
			
          \((TłAutoK[[i, 
                PrevMaxLoc[[2
                  ]]]])\)*
            \((PrevTable
                [[PrevMaxLoc[[2
                  ]]]])\); 
	
      AppendTo[sfa, sfal]; 
\ \ 
      PrevTable = TłAutoK[[i]]; 
      
	PrevMaxLoc = MaxLoc]; 
sfa *= 4; \)

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

Parametr segmentacji

SegMeas = (Sqrt[(sfa^2) + (sft^2)]);

Miara segmentacji

ListPlot[SegMeas, PlotRange -> {0, 2},  PlotJoined -> True, 
      PlotStyle -> RGBColor[1, 0, 0], Frame -> True];

Próg segmentacji

SegTable = 
      Table[SegMeas[[i]] < 1, {i, 2, 
          Length[SegMeas]}];

Tabela punktów segmentacji

SegPoints = 
      Table[If[SegTable[[i]] && 
            ! SegTable[[i - 1]], 1, 
          0], {i, 2, Length[SegTable]}];

\(Gap = 10; 

    For[i = 1, i <= Length[SegPoints] - Gap, \(i++\), 
	
      If[SegPoints[[i]] == 1, 
		
        For[j = 1, j <= Gap - 1, \(j++\), 
			
          SegPoints[[i + j]] = 0]]]; 
    \)

Liczba punktów podziału

SegN = 
      1 + (∑_(i = 1)^(Length[SegPoints]))SegPoints
            [[i]]

Analiza sygnału w segmentach

Obliczenie redniej ważonej funkcji autokorelacji w segmentach

\(j = 0; 
prevj = 0; 
MeanWall = {}; 

    MeanWAK = 
      Table[0, {
          Length[TłAutoK[[1
              ]]]}]; 
Wy = {}; 

    For[i = 1, i <= SegN, \(i++\), 
	AutoKTable = {}; 
	\(j++\); 
	
      While[SegPoints[[j]] != 1, 
		
        AppendTo[AutoKTable, 
          TłAutoK[[j]] 
            SegMeas[[j]]]; \(j++\)]; 
	
      For[k = 1, 
        k <= Length[
            TłAutoK[[1]]], 
        \(k++\), 
		
        MeanWAK[[k]] = 
          Mean[Table[
              AutoKTable[[l, k]], {l, 
                j - prevj - 1}]]]; 
	prevj = j; 
	
      AppendTo[MeanWall, MeanWAK\/Max[MeanWAK]]]; \)

Detekcja maksimów funkcji autokorelacji

MaksSegAK = MaksAutoK[MeanWall];

Interfejs predyktora

Obliczenie długoci segmentów

\(SegL = {}; 
prev = 0; 

    For[i = 1, i <= Length[SegPoints], \(i++\), 
	
      If[SegPoints[[i]] == 1, 
		
        AppendTo[SegL, i - prev]; 
		prev = i]]; 

    AppendTo[SegL, Length[SegPoints] - prev]; \)

Wyznaczenie zakresów poszczególnych pasm półtonowych

Pasma = GetBands[444, CzMin, CzMax];

Wyznaczenie maksimów funkcji autokorelacji w pasmach dla poszczególnych segmentów

MaxWPasmach = GetMaksBands[MeanWall, Pasma];

ListPlot3D[MaxWPasmach];

Zapis danych do pliku

\(OutName = "autoc_"<>ToString[Tłumik]<>".data"; 

    OutStream = OpenWrite[OutName]; 

    WriteString[OutStream, 
      "Autocorrelation data
"<>ToString[Length[MaxWPasmach]]<>" "<>
        ToString[
          Length[MaxWPasmach[[1]]]]<>
        "
"]; 

    For[i = 1, i <= Length[MaxWPasmach], \(i++\), 
	
      For[j = 1, 
        j <= Length[MaxWPasmach[[i]]], 
        \(j++\), 
		
        WriteString[OutStream, 
          N[MaxWPasmach[[i, j]]]]; 
		
        If[j < Length[MaxWPasmach[[i]]], 
          
			WriteString[OutStream, " "], 
			
          WriteString[OutStream, 
            " "<>
              ToString[SegL[[i]]]<>
              "
"]]]]; 
Close[OutStream]; \)

Odczyt danych z pliku

\(InName = "D:\source\mtrack\output"; 

    Res = ReadList[InName, \ Number]; \)

CzSeg = 
      ((Pasma[[Res]] + 
          Pasma[[Res + 1]])/2);

MaksSegAK = Round[(CzP/CzSeg)];

Analiza wyników

Wywietlenie wykresów

ListPlot[
      Table[(CzP/MaksSegAK)[[x]], {x, 
          Length[MaksSegAK]}], PlotRange -> {{0, SegN + 1}, All}, 
      PlotStyle -> RGBColor[0, 0, 1], Frame -> True, Axes -> False];

Analiza pojedynczych ramek sygnału

Znalezienie maksimum lokalnego funkcji autokorelacji położonego najbliżej (Subscript[f, 0]) wyznaczonej w segmencie

\(TabFreq = {}; 
TabVal = {}; 
prevloc = 0; 
k = 1; 

    For[i = 1, i <= Length[TłAutoK], \(i++\), 
	
      If[SegPoints[[i]] == 1, \(k++\)]; 
      
	AKVal = 
        Max[TłAutoK[[i]]]; 
	
      loc = \(initloc = 
          MaksSegAK[[k]]\); 
	
      While[TłAutoK[[i, loc]] < 
            TłAutoK[[i, 
              loc + 1]] || 
          TłAutoK[[i, loc]] < 
            TłAutoK[[i, 
              loc - 1]], 
		
        If[TłAutoK[[i, 
              loc + 1]] > 
            TłAutoK[[i, 
              loc - 1]], \(loc++\), \(loc--\)]]; 
	
      ToF = {TłAutoK[[i, 
            loc - 1]], 
          TłAutoK[[i, loc]], 
          TłAutoK[[i, 
            loc + 1]]}; 
	
      If[SegMeas[[i]] < 1 && 
          \((loc\/initloc > \@2\%12 || loc\/initloc < \@\(2\^\(-1\)\)\%12)\), 
        loc = prevloc]; 
	prevloc = loc; 
	
      poly = \(-Fit[ToF, {1, x, x\^2}, x]\); 
	
      korekta = 
        x /. \(FindMinimum[poly, {x, 1}]\)[[2
            ]]; 
	acloc = loc + korekta - 2; 
	
      AppendTo[TabFreq, CzP\/acloc]; 
]; \)

ListPlot[TabFreq, PlotJoined -> False, 
      PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.005]}, Frame -> True, 
      Axes -> False];

Wykres porównawczy z danymi MIDI

Przygotowanie wykresu

\(Wykresy = {}; 

    For[i = 1, i <= Length[TabFreq], \(i++\), 
	
      AppendTo[Wykresy, 
        Graphics[{RGBColor[0, 0, 1], PointSize[0.004], 
            Point[{\(i*Skok\)\/CzP, 
                TabFreq[[i]]}]}]]]; \)

Show[{Graph, Wykresy}, Frame -> True, Axes -> False];

Analiza błędów

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

\(\(Błędzik[\([2]\)] = 0; \)\)

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

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