Detekcja częstotliwoci podstawowej przy pomocy analizy obwiedni sygnału

Š 1999, 2000, Marek Szczerba

Inicjalizacja

Wczytanie bibliotek

Biblioteka obsługi dwięku

<< Miscellaneous`Audio`

Biblioteki statystyczne

<< Statistics`DescriptiveStatistics`

Akwizycja sygnału

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

snd = ReadSoundfile["c3forts.wav", PrintHeader -> True];

Parametry definiowane przez użytkownika

Parametr tłumienia

P=(10/1000);

Częstotliwoć próbkowania

CzP=22050;

Granice pasma analizy

CzMin=100;
CzMax=2000;

Wysokoć dwięku w skali muzycznej (różnica od (a^1))

Wys=15;

Parametry obliczane automatycznie

Konwersja pasma analizy

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

Rzeczywista częstotliwoć podstawowa dwięku

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

Definicje funkcji przetwarzania sygnału

Detekcja maksimów i pominięcie pozostałych próbek

\(\(Maks = 
      Function[dane, 
		Wynik = Table[0, {fi, Length[dane]}]; 
		
        For[fi = 2, fi < Length[dane], \(fi++\), 
			
          \(If[Abs[dane[[fi]]] > 
                Abs[dane[[fi - 1]]] && 
              Abs[dane[[fi]]] > 
                Abs[dane[[fi + 1]]], 

            					
            \(Wynik[[fi]] = 
              dane[[fi]]; \)]; \)]; 

        		Wynik]; \)\)

Generowanie obwiedni

Obwiednia jednostronna (dodatnia)

\(\(\(Obw1Str = 
      Function[{dane, wsp}, 
		PlusMaks = 0; 
		
        Env = Table[0, {fn, Length[dane]}]; 
		
        For[fn = 1, fn <= Length[dane], \(fn++\), 
			
          If[dane[[fn]] > PlusMaks && 
              dane[[fn]] > 0, 
				
            PlusMaks = dane[[fn]], 

            				PlusMaks = PlusMaks*Exp[\(-wsp\)]]; 
		
          Env[[fn]] = N[PlusMaks]; ]; 

        		Env]; \)			\)\)

Obwiednia dwustronna

\(\(Obw2Str = 
      Function[{dane, wsp}, 
		MinusMaks = \(PlusMaks = 0\); 
		
        LowEnv = \(HighEnv = Table[0, {fn, Length[dane]}]\); 
		
        For[fn = 1, fn <= Length[dane], \(fn++\), 
			
          If[dane[[fn]] > PlusMaks && 
              dane[[fn]] > 0, 
				
            \(PlusMaks = dane[[fn]]; \), 
            
				PlusMaks = PlusMaks*Exp[\(-wsp\)]]; 
			
          If[PlusMaks < 0, PlusMaks = 0]; 
			
          If[dane[[fn]] < MinusMaks && 
              dane[[fn]] < 0, 
				
            \(MinusMaks = dane[[fn]]; 
            \), 
				MinusMaks = MinusMaks*Exp[\(-wsp\)]]; 
			
          If[MinusMaks > 0, MinusMaks = 0]; 
			
          LowEnv[[fn]] = N[MinusMaks]; 

          			HighEnv[[fn]] = 
            N[PlusMaks]; ]; 
	Wynik = {LowEnv, HighEnv}; 
	Wynik]; \)\)

Redukcja obwiedni

\(\(RedObw = 
      Function[dane, 
		Wynik = Table[0, {fn, Length[dane]}]; 
		
        Wynik[[1]] = 
          dane[[1]]; 
		
        For[fn = 2, fn <= Length[dane], \(fn++\), 
			
          If[dane[[fn]] > 
              dane[[fn - 1]], 
				
            Wynik[[fn]] = 
              dane[[fn]] - 
                dane[[fn - 1]]]]; 
	
        Wynik]; \)\)

Detekcja poszczególnych okresów sygnału

\(\(Okr = 
      Function[dane, 
		Wynik = {}; 
		Poprz = 0; 
		Amp = 0; 
		
        For[fn = 1, fn < Length[dane], \(fn++\), 
			
          If[Sign[dane[[fn]]] > 0, 

            				Amp = dane[[fn]]; 

            				If[fn - Poprz >= OkrMin, 
					
              If[fn - Poprz > OkrMax, Amp = 0]; 
					
              AppendTo[Wynik, {fn, Amp}]; 
					Poprz = fn; ]]]; 
		
        Wynik]; \)\)

Obliczenie poprawnoci detekcji

\(\(BłšdProc = 
      Function[dane, 
		GrDolna = CzP\/\(RzCzPodst \@\( 2\^\(-1\)\)\%12\); 
        
		GrGórna = CzP\/\(RzCzPodst \@ 2\%12\); 
		Wynik = 0; 
        
		Poprz = 0; 
		
        For[fn = 1, fn < Length[dane], \(fn++\), 
			
          If[dane[[fn, 1]] - Poprz < 
                GrDolna && 
              dane[[fn, 1]] - Poprz > 
                GrGórna, \(Wynik++\)]; 
			
          Poprz = dane[[fn, 1]]; ]; 

        		If[Length[dane] == 0, 
			Wynik = 0, 
			
          \(Wynik = N[100  Wynik\/Length[dane]]; \)]; 
		Wynik]; \)\)

Obliczenie błędu ledzenia częstotliwoci podstawowej (w półtonach na okres)

\(\(FBłędu = 
      Function[dane, 
			Wynik = {}; 
			Poprz = 0; 
			
        For[fn = 1, fn < Length[dane], \(fn++\), 
				
          If[dane[[fn, 2]] > 0, 

            					
            AppendTo[Wynik, 
              CzP\/\(dane[[fn, 1]] - 
                    Poprz\) - RzCzPodst]]; 
				
          Poprz = dane[[fn, 1]]; ]; 

        			Wynik]; \)\)

Obliczenia

Detekcja maksimów i pominięcie próbek

Maksima = Maks[snd];

Obwiednia z tłumieniem

Obwiednia = Obw1Str[Maksima, P];

Detekcja maksimów obwiedni

Redukcja = RedObw[Obwiednia];

Obliczenie obwiedni maksimów obwiedni

Obwiednia2 = Obw1Str[Redukcja, P];

Detekcja maksimów obwiedni

Redukcja2 = RedObw[Obwiednia2];

Estymacja poszczególnych okresów sygnału

Okresy = Okr[Redukcja2];

Obliczenie błędu detekcji

\(\(Błšd = FBłędu[Okresy]; \)\)

Wykres błędu

ListPlot[Błšd, PlotStyle -> RGBColor[0, 0, 1], 
      PlotJoined -> True, PlotRange -> All];

rednia wartoć błędu

N[Mean[Abs[Błšd]]]

Poprawnoć detekcji

Epr = BłšdProc[Okresy]

Przygotowanie tabeli zawierajšcej wykresy częstotliwoci podstawowej poszczególnych okresów sygnału

\(Wykresy = {}; 

    For[i = 1, i < Length[Okresy], \(i++\), 
	
      If[Okresy[[i, 2]] > 0, 
		
        AppendTo[Wykresy, 
          Graphics[{RGBColor[0, 0, 1], PointSize[0.004], 
              Point[{Okresy[[i + 1, 1
                      ]]\/CzP, 
                  CzP\/\(Okresy[[i + 1, 1
                        ]] - 
                      Okresy[[i, 1
                        ]]\)}]}]]]]; \)

Wywietlenie wykresów

Show[Wykresy, Frame -> True, PlotRange -> All];

Rekonstrukcja sygnału

\(Wy = {}; 
Poprz = 0; 

    For[i = 1, i <= Length[Okresy], \(i++\), 
	
      A = Sign[Okresy[[i, 2]]]; 
	
      AppendTo[Wy, 
        Table[A\ Sin[j], {j, 
            N[\(2  π\)\/\(Okresy[[i, 1
                    ]] - Poprz\)], 2  π, 
            N[\(2  π\)\/\(Okresy[[i, 1
                    ]] - Poprz\)]}]]; 
	
      Poprz = Okresy[[i, 1]]; ]; 
	
    Wy = Flatten[Wy]; \)

ListPlay[Wy, SampleRate -> CzP, PlayRange -> All];

Optymalizacja parametru tłumienia obwiedni

\(\(P = 1\/30; 
R = 2; 
Niegotowe = True; 

    While[Niegotowe, 
	Parametry = {P - P\/R, P, P + P\/R}; 
	
      OK = {0, 0, 0}; 
	
      For[i = 1, i <= Length[Parametry], \(i++\), 
		
        Obwiednia = 
          Obw1Str[Maksima, 
            Parametry[[i]]]; 
		
        Redukcja = RedObw[Obwiednia]; 
		
        Obwiednia2 = 
          Obw1Str[Redukcja, 
            Parametry[[i]]]; 
		
        Redukcja2 = RedObw[Obwiednia2]; 
		Okresy = Okr[Redukcja2]; 
		
        OK[[i]] = 
          BłšdProc[Okresy]; 
	]; 
	
      Print["P=", N[P], "   R=", R, "   OK=", OK]; 
	
      If[Max[OK] == OK[[2]], R = 2  R, 

        		If[OK[[1]] > 
            OK[[3]], 
          P = Parametry[[1]], 
          P = Parametry[[3]]]]; 
	
      If[R > 100, Niegotowe = False]; ]; 			\)\)