Histogram Schroedera

Š Marek Szczerba, 2000

Inicjalizacja

Wczytanie bibliotek

<< Miscellaneous`Audio`<< Graphics`Graphics`<< CCRMA`SCMTheory`<< CCRMA`Windows`<< Statistics`MultiDescriptiveStatistics`;

Akwizycja sygnału

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

snd = ReadSoundfile["emfmobs.wav"];

Parametry definiowane przez użytkownika

Częstotliwoć próbkowania

CzP=22050;

Rozmiar okna czasowego

OknoN=512;

Zakładka

Zakładka=1/2;

Liczba harmonicznych do analizy

MaksHarm=7;

Współczynnik zeropaddingu

wZP=3;

Zakres częstotliwoci podstawowej

CzMin=100;
CzMaks=2000;

Próg szumowy (w dB)

PrSz=-25;

Próg amplitudy maksimów (w dB względem amplitudy max. pršżka w widmie)

PrAmp=-20;

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

        	]; \)\)

Detekcja maksimów widma

\(\(MaksWidm = 
      Function[{dane, wsp, próg}, 
		fMaksy = {}; 
		Wy = {}; 

        		fmax = 0; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			
          If[dane[[fi]] > fmax, 
            fmax = dane[[fi]]]]; 
		
        fminamp = fmax\ 2\^\(próg\/3\); 
		
        For[fi = 2, fi < Length[dane], \(fi++\), 
			
          If[dane[[fi]] > 
                dane[[fi - 1]] && 
              dane[[fi]] > 
                dane[[fi + 1]] && 
              dane[[fi]] > fminamp, 

            				AppendTo[
              fMaksy, {dane[[fi]], 
                fi}]]]; 
		fMaksy = Sort[fMaksy, Greater]; 
		
        If[Length[fMaksy] > wsp, 
			
          For[fi = 1, fi <= wsp, \(fi++\), 
            AppendTo[Wy, fMaksy[[fi]]]], 
          Wy = fMaksy]; 
		
        Wy = Table[{Wy[[fi, 2]], 
              Wy[[fi, 1]]}, {fi, 
              Length[Wy]}]; 
	Sort[Wy]]; \)\)

Histogram Shroedera

\(\(Histogram = 
      Function[{dane}, 
		
        Loc = Table[
            dane[[x, 1]], {x, 
              Length[dane]}]; 
		
        MaxRel = 
          Ceiling[Loc[[Length[Loc]]]\/Loc
                [[1]]]; 
		
        Wy = Table[0, {x, Length[Loc]}]; 
		
        For[fi = 1, \ fi < Length[Loc], \(fi++\), 
			
          For[fj = 2, \ fj <= MaxRel, \ \(fj++\), 
			
            For[fk = fi + 1, fk <= Length[Loc], \(fk++\), 
				
              If[Abs[Loc[[fk]]\/fj - 
                        Loc[[fi]]] < 
                    1 || Abs[
                      \(Loc[[fk]] - 1
                            \)\/fj - 
                        Loc[[fi]]] < 
                    1 || Abs[
                      \(Loc[[fk]] + 1
                            \)\/fj - 
                        Loc[[fi]]] < 1, 
                \(Wy[[fi]]++\)]]]]; 	Wy
        
		]; \)\)

Detekcja maksimów

\(\(FMaks = 
      Function[{dane, fwsp}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			fmw = 0; fml = 0; 

          			For[fj = 1, 
            fj <= Length[dane[[fi]]], 
            \(fj++\), 
				
            If[dane[[fi, fj]] > fmw, 

              					
              fmw = dane[[fi, fj]]; 

              					fml = fj; ]]; 
				
          AppendTo[Wy, {fmw, fml\ fwsp}]]; 
		Wy]; \)\)

Detekcja maksimum histogramu i zestawienie z widmem

\(\(Dop = 
      Function[{fwidmo, fhist}, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[fhist], \(fi++\), 
			fmax = {0, 0}; 

          			For[fj = 1, 
            fj <= Length[fhist[[fi]]], 
            \(fj++\), 
				
            If[fhist[[fi, fj]] > 
                fmax[[2]], 
              fmax = {fj, 
                  fhist[[fi, fj]]}]]; 

          			If[fmax[[2]] > 0, 
            AppendTo[Wy, 
              fwidmo[[fi, 
                fmax[[1]], 1
                ]]], AppendTo[Wy, 0]]]; 
	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)\)]; ]; \)

Detekcja maksimów widma

\(MaksimaWidma = {}; 

    For[i = 1, i <= Length[Widma], \(i++\), 
      AppendTo[MaksimaWidma, 
        Sort[MaksWidm[Widma[[i]], 
            MaksHarm, PrAmp]]]]; \)

Obliczenie histogramu Shroedera

\(Hist = {}; 

    For[i = 1, i <= Length[MaksimaWidma], \(i++\), 
      AppendTo[Hist, 
        Histogram[MaksimaWidma[[i]]]]]; 
    \)

Detekcja maksimów histogramu i dopasowanie częstotliwoci

Fs = Dop[MaksimaWidma, Hist] (CzP/((1 + wZP) OknoN));

Konwersja czasu

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

Wykres błędu

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

Show[a, Graph];

Wartoć rednia miary błędu

Dobre = BłšdProc[Fs]

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

Optymalizacja parametru MaksHarm

\(P = 2; 
NWynik = 0; 
niegotowe = True; 

    While[niegotowe, 
	MaksimaWidma = {}; 
	
      For[i = 1, i <= Length[Widma], \(i++\), 
        AppendTo[MaksimaWidma, 
          Sort[MaksWidm[Widma[[i]], P, 
              PrAmp]]]]; 
	Hist = {}; 
	
      For[i = 1, i <= Length[MaksimaWidma], \(i++\), 
        AppendTo[Hist, 
          Histogram[
            MaksimaWidma[[i]]]]]; 
	
      Fs = Dop[MaksimaWidma, Hist] CzP\/\(\((1 + wZP)\) OknoN\); 
	
      Dobre = BłšdProc[Fs]; 
	
      Print["P=", P, "   OK=", Dobre]; 
	\(P++\); 
	
      If[Dobre >= NWynik, NWynik = Dobre, If[P > 15, niegotowe = False]]; ]; 
    \)