Metoda grzebieniowa

Š Marek Szczerba, 2000

Inicjalizacja

Wczytanie bibliotek

<< Miscellaneous`Audio`<< Graphics`Graphics`\(<< Graphics`Graphics3D`
<< CCRMA`SCMTheory`\)<< CCRMA`Windows`<< Statistics`MultiDescriptiveStatistics`;

Akwizycja sygnału

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

    snd = ReadSoundfile["mflets.wav"]; \)

Parametry okrelane przez użytkownika

Częstotliwoć próbkowania

CzP=22050;

Rozmiar okna czasowego

OknoN=512;

Zakładka

Zakładka=1/2;

Współczynnik zeropaddingu

wZP=3;

Próg szumowy (w dB)

PrSz=-∞;

Zakres częstotliwoci podstawowej

CzMin=100;
CzMax=2000;

Wartoć tłumienia wyższych częstotliwoci grzebienia

Tłumik=2;

Detekcja maksymalnej wartoci bezwzględnej próbki w sygnale

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

Zakres częstotliwoci podstawowej (w próbkach)

\(FMin = Floor[CzMin \((\((1 + wZP)\)\ OknoN)\)/CzP]; 

    FMax = Ceiling[CzMax \((\((1 + wZP)\)\ OknoN)\)/CzP]; \)

Definicje funkcji przetwarzania sygnału

Definicja okna czasowego

\(\(FOkna = Function[dane, 
		dane*HammingW[Length[dane]]]; \)\)

Nadpróbkowanie z interpolacjš wielomianem 2-go stopnia

\(\(NadPr = 
      Function[{dane, wsp}, 
		ovr = Table[0, {x, Length[dane] wsp}]; 

        		For[fi = 2, fi < Length[dane], \(fi++\), 
			
          fpróbki = 
            Table[dane[[fi + fj]], {fj, 
                \(-1\), 1}]; 
			
          fpoly = Fit[fpróbki, {1, x, x\^2}, x]; 
			
          For[fj = 0, fj < wsp, \(fj++\), 
				
            ovr[[fi\ wsp + fj]] = 
              fpoly /. x -> \((1 + fj/wsp)\)]; ]; 
			ovr]; \)\)

Zeropadding

\(\(ZP = 
      Function[{dane, wZP}, 
		Wy = dane; 
		
        zera = Table[0, {x, wZP\ Length[dane]}]; 
		AppendTo[Wy, zera]; 

        		Flatten[Wy]]; \)\)

Obliczenie widm

\(\(OblWidma = 
      Function[dane, 
		TempWy = Abs[Fourier[dane]]; 
		
        Wy = Table[
            TempWy[[fx]], {fx, 
              Length[TempWy]\/2}]; 
		
        For[fi = 1, fi <= FMin, \(fi++\), 
          \(Wy[[fi]] = 0; \)]; 
		Wy

        	]; \)\)

Definicja funkcji grzebieniowych bez tłumienia wyższych składowych

\(\(Grz1 = 
      Function[fl, 
		Wy = {}; 
		
        For[fi = FMin, fi <= FMax, \(fi++\), 
			
          fComb = Table[0, {x, fl}]; 
			
          For[fj = fi, fj <= fl, fj += fi, 
				
            fComb[[fj]] = 1]; 
			
          AppendTo[Wy, fComb]]; 
	Wy]; \)\)

Definicja funkcji estymacji harmonicznej z tłumieniem wyższych częstotliwoci

\(\(Grz2 = 
      Function[{fl, wsp}, 
		Wy = {}; 
		
        For[fi = FMin, fi <= FMax, \(fi++\), 
			
          fComb = Table[0, {x, fl}]; 
			
          For[fj = 1, fj*fi < fl, \(fj++\), 
				
            fComb[[fj*fi]] = 
              N[fj\^\(-\(1\/wsp\)\)]]; 
			AppendTo[Wy, fComb]]; 
			
        Wy
	]; \)\)

Definicja funkcji estymacji harmonicznej

\(\(HEF = 
      Function[{X, G}, 
		Wy = Table[0, {FMax}]; 
		
        For[fi = FMin, fi <= FMax, \(fi++\), 
			
          Wy[[fi]] = 
            ∑\+\(fk = 1\)\%\(Length[X]\/fi\)X
                  [[fi\ fk]] 
                G[[fi - FMin + 1, 
                  fi\ fk]]]; 
				Wy]; \)\)

Definicja funkcji detekcji maksimów

\(\(Maksima = 
      Function[dane, 
		Wy = Table[0, {Length[dane]}]; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			fMaks = {0, 0}; 

          			For[fj = 1, 
            fj <= Length[dane[[fi]]], 
            \(fj++\), 
				
            If[dane[[fi, fj]] > 
                fMaks[[1]], 
              fMaks[[1]] = 
                dane[[fi, fj]]; 
              fMaks[[2]] = fj]]; 
			
          Wy[[fi]] = 
            fMaks[[2]]]; 
		Wy]; 
    \)\)

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)\) OknoN\ \((1 - Zakładka)\)\)\/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 widma

\(i = 0; 
Widma = {}; 

    While[i <= Length[snd] - OknoN, 
	
      Ramka = Table[
          snd[[x]], {x, i + 1, 
            i + OknoN}]; 
	
      If[20  Log[10, Max[Abs[Ramka]]\/MaksPróbka] > PrSz, 
		
        Okno = FOkna[Ramka]; 
		ZeroP = ZP[Okno, wZP]; 
		
        AppendTo[Widma, OblWidma[ZeroP]]]; 

      i += Floor[OknoN \((1 - Zakładka)\)]; ]; \)

Wygenerowanie funkcji grzebieniowych

GrzebTł = 
      Grz2[Length[Widma[[1]]], 
        Tłumik];

Obliczenie funkcji estymacji harmonicznej

\(FEH = {}; 

    For[i = 1, i <= Length[Widma], \(i++\), 
	
      AppendTo[FEH, 
        HEF[Widma[[i]], \ 
          GrzebTł]]]; \)

Detekcja maksimów i przeliczenie na Hz

FS = Maksima[FEH] (CzP/(OknoN (1 + wZP)));

Konwersja czasu

InS = 
      Table[{(((i - 1) OknoN (1 - Zakładka))/CzP), 
          FS[[i]]}, {i, Length[FS]}];

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

Show[a, Graph];

Obliczenie błędu

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

Wykres błędu

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

Sprawdzenie funkcjonowania metody dla różnych wartoci parametru s

\(Dobre = {}; 
Błędy = {}; 

    For[P = 1, P <= 10, \(P++\), 
	Print[P]; 
	
      GrzebTł = 
        Grz2[Length[Widma[[1]]], P]; 
	
      FEH = {}; 
	
      For[j = 1, j <= Length[Widma], \(j++\), 
		
        AppendTo[FEH, 
          HEF[Widma[[j]], \ 
            GrzebTł]]]; 
	
      FS = Maksima[FEH] CzP\/\(OknoN \((1 + wZP)\)\); 
	
      AppendTo[Dobre, BłšdProc[FS]]; 
	
      AppendTo[Błędy, Mean[Abs[Błędzik]]]; 
]; 
    \)

\(\({Dobre, Błędy} // N\) // MatrixForm\)