Detektor Beauchampa

Š Marek Szczerba, 2000

Inicjalizacja

Wczytanie bibliotek

\(<< Miscellaneous`Audio`
<< CCRMA`SCMTheory`
<< CCRMA`Windows`

      \(<< Statistics`MultiDescriptiveStatistics`; \)\)

Akwizycja sygnału

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

    snd = ReadSoundfile["c2orgs.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;

Zakres częstotliwoci

CzMin=100;
CzMaks=2000;

Maksymalna liczba analizowanych harmonicznych

MaksHarm=10;

Próg szumowy (w dB)

PrSz=-25;

Próg amplitudowy detekcji maksimów widma

PrAmp=-30;

Waga błędu (Subscript[e, 2])

W=3;

Zakres częstotliwoci podstawowej

CzMin=100;
CzMax=2000;

Detekcja maksymalnej wartoci bezwzględnej próbki w sygnale

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

Rzeczywista częstotliwoć podstawowa

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

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]]]; \)\)

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

        	]; \)\)

Detekcja maksimów widma i interpolacja wielomianem 2-stopnia

\(\(DetIntMaks = 
      Function[{dane, wsp, licz}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			fmaks = 0; 
			
          For[fj = FMin, 
            fj <= Length[dane[[fi]]], 
            \(fj++\), 
				
            If[dane[[fi, fj]] > fmaks, 
              fmaks = dane[[fi, fj]]]]; 
          
			fFr = {}; 
			
          For[fj = FMin, fj <= FMax, \(fj++\), 
				
            If[dane[[fi, fj]] > 
                  dane[[fi, fj + 1]] && 
                dane[[fi, fj]] > 
                  dane[[fi, fj - 1]] && 
                20  Log[10, 
                      dane[[fi, fj
                          ]]\/fmaks] > wsp, 
					
              fpoly = Fit[{
                    dane[[fi, fj - 1]], 
                    dane[[fi, fj]], 
                    dane[[fi, 
                      fj + 1]]}, {1, x, x\^2}, x]; 

              					
              AppendTo[
                fFr, {dane[[fi, fj
                      ]]\/fmaks, 
                  fj - 2 + x /. 
                    \(FindMinimum[\(-fpoly\), {x, 0}]\)[[2
                      ]]}]]]; 
			
          If[licz < Length[fFr], 
				fFr = Sort[fFr, Greater]; 

            				fFr = 
              Table[fFr[[fj]], {fj, 
                  licz}]]; 
			
          fFr = Table[{fFr[[fj, 2]], 
                fFr[[fj, 1]]}, {fj, 
                Length[fFr]}]; 
			AppendTo[Wy, fFr]]; 
		Wy]; \)\)

Obliczenie błędu dopasowania (Subscript[e, 1])

\(\(OblE1 = 
      Function[{dane, wsp}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			eFr = {}; 
			
          For[fh = FMin, fh <= FMax, fh = N[fh \@ 2\%12], 
				fe = 0; 

            				fdist = ∞; 
				
            For[fj = 1, fj <= wsp, \(fj++\), 
					
              For[fk = 1, 
                fk <= Length[
                    dane[[fi]]], 
                \(fk++\), 
						
                If[Abs[dane[[fi, fk, 1
                          ]] - fh*fj] < fdist, 
                  fdist = 
                    Abs[dane[[fi, fk, 1
                          ]] - fh*fj]]]; 
					
              fe += fdist]; 
				AppendTo[eFr, fe]]; 
			
          AppendTo[Wy, eFr]]; 
		Wy]; \)\)

Obliczenie błędu dopasowania (Subscript[e, 2])

\(\(OblE2 = 
      Function[{dane, wsp}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			eFr = {}; 
			
          For[fh = FMin, fh <= FMax, fh = N[fh \@ 2\%12], 
				fe = 0; 

            				fdist = ∞; 
				
            For[fk = 1, 
              fk <= Length[dane[[fi]]], 
              \(fk++\), 
					fhfj = fh; 
					
              For[fj = 1, fj <= wsp, \(fj++\), 
						fhfj = fh\ fj; 

                						
                If[Abs[dane[[fi, fk, 1
                          ]] - fhfj] < fdist, 
                  fdist = 
                    Abs[dane[[fi, fk, 1
                          ]] - fhfj]]]; 
					
              fe += fdist]; 
				AppendTo[eFr, fe]]; 
			
          AppendTo[Wy, eFr]]; 
		Wy]; \)\)

Funkcja detekcji minimów i konwersji do Hz

\(\(DetMinK = 
      Function[dane, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			
          fmin = {∞, 0}; 
			
          For[fj = 1, 
            fj <= Length[dane[[fi]]], 
            \(fj++\), 
				
            If[dane[[fi, fj]] < 
                fmin[[1]], 
              fmin = {dane[[fi, fj]], 
                  fj}]]; 
			
          AppendTo[Wy, 
            N[FMin \( CzP\/\(OknoN \((1 + wZP)\)\)\) 
                \@\(2\^fmin[[2
                        ]]\)\%12]]]; 
		Wy]; \)\)

Definicja funkcji błędu

\(\(FBłędu = 
      Function[{dane, rzcz}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			
          If[dane[[fi]] > 0, 
            AppendTo[Wy, 
              dane[[fi]] - rzcz]]]; 

        		Wy]; \)\)

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)\)]; ]; \)

Detekcja maksimów widma

Maks = DetIntMaks[Widma, PrAmp, MaksHarm];

Obliczenie błędu dopasowania (Subscript[e, 1])

e1 = OblE1[Maks, 10];

e2 = OblE2[Maks, 10];

ListPlot3D[e1 + e2, PlotRange -> {All, All, {0, 200}}];

Obliczenie błędu sumarycznego E

ES = e1 + W e2;

Detekcja minimów i konwersja na Hz

FS = DetMinK[ES];

Wykres wyniku detekcji

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

Obliczenie błędu detekcji

\(\(Błšd = FBłędu[FS, RzCzPodst]; \)\)

Wykres błędu

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

Wartoć rednia miary błędu

Mean[Abs[Błšd]]

Wartoć miary poprawnoci detekcji częstotliwoci podstawowej

N[DFpr[FS, RzCzPodst]]

Badanie dla różnych wartoci wagi W

\(Błędy = {}; 

    For[i = \(-10\), i <= 10, \(i++\), 
	ES = e1 + 2^i\ e2; 
	
      FS = DetMinK[ES]; 
	
      Błšd = FBłędu[FS, RzCzPodst]; 
	
      AppendTo[Błędy, {i, Mean[Abs[Błšd]], 
          N[DFpr[FS, RzCzPodst]]}]]; \)

\(Błędy\)