Model Meddisa - Hewitta

Š Marek Szczerba, 2000

Inicjalizacja

Wczytanie bibliotek

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

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

Akwizycja sygnału

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

    snd = ReadSoundfile["emfmoboe2.wav", PrintHeader -> True]; \)

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=0;

Zakres częstotliwoci

CzMin=100;
CzMaks=2000;

Próg szumowy (w dB)

PrSz=-∞;

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[CzMaks \((\((1 + wZP)\)\ OknoN)\)/CzP]; \)

Pasma częstotliwoci

\(Pasma = {}; 
fx = Floor[\(50 \((1 + wZP)\) OknoN\)\/CzP]; 

    fmax = Ceiling[\(6000 \((1 + wZP)\) OknoN\)\/CzP]; 

    While[fx < fmax, 
	fy = fx \@\( 2\^8\)\%12; 
	
      If[fy > fmax, fy = fmax]; 
	AppendTo[Pasma, {fx, Ceiling[fy]}]; 
	
      If[fx == Floor[fy], \(fx++\), fx = Floor[fy]]]; \)

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

Funkcja okna czasowego i zeropaddingu dla pasma

\(\(ZeroPHamm = 
      Function[danew, 
		Wyj = {}; 
		
        For[fi = 1, fi <= Length[danew], \(fi++\), 
			
          framka = FOkna[danew[[fi]]]; 

          			AppendTo[Wyj, ZP[framka, 3]]]; 
		Wyj]; \)\)

Obliczenie widm

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

        	]; \)\)

Funkcja podziału widma na pasma

\(\(Podział = 
      Function[{dane, pasma}, 
		
        Wy = Table[Table[0, {Length[dane]}], {Length[pasma]}]; 
		
        For[fj = 1, fj <= Length[pasma], \(fj++\), 
			
          For[fk = pasma[[fj, 1]], 
            fk <= pasma[[fj, 2]], 
            \(fk++\), 
            Wy[[fj, fk]] = 
              dane[[fk]]]]; 
		Wy]; 
    \)\)

Rekonstrukcja sygnału w pasmach

\(\(RekPasm = 
      Function[{dane}, 
		Wy = {}; 
		
        For[fj = 1, fj <= Length[dane], \(fj++\), 
			
          fwidmo = 
            Table[0, {
                Length[dane[[fj]]] 2}]; 
          
				
          For[fk = 1, 
            fk <= Length[dane[[fj]]], 
            \(fk++\), 
					
            fwidmo[[fk]] = 
              dane[[fj, fk]]; 

            					
            fwidmo[[Length[fwidmo] + 1 - fk
                  ]] = 
              dane[[fj, fk]]]; 
				
          fsygn = Re[InverseFourier[fwidmo]]; 
				AppendTo[Wy, fsygn]]; 
        
		Wy]; \)\)

Prostownik pół fali

\(\(Prostownik = 
      Function[dane, 
		Wy = dane; 
		
        For[fj = 1, fj <= Length[Wy], \(fj++\), 
			
          For[fk = 1, 
            fk < Length[Wy[[fj]]], 
            \(fk++\), 
				
            If[dane[[fj, fk]] < 0, 
              Wy[[fj, fk]] = 0]]]; 
		
        Wy]; \)\)

Funkcja obwiedni amplitudy

\(\(FOAM = 
      Function[{dane, pasma}, 
		Wy = {}; 
		
        For[fj = 1, fj <= Length[pasma], \(fj++\), 
			
          fwidmo = Fourier[dane[[fj]]]; 
          
			For[
            fk = pasma[[fj, 1]], 
            fk < Length[fwidmo] - 
                pasma[[fj, 1]] + 1, 
            \(fk++\), fwidmo[[fk]] = 0]; 
          
			fsygn = Re[InverseFourier[fwidmo]]; 
			
          AppendTo[Wy, fsygn]]; 
		Wy]; \)\)

Estymacja funkcji autokorelacji

\(\(ACF = 
      Function[dane, 
		Wy = {}; 
		
        For[fi = 1, fi <= Length[dane], \(fi++\), 
			
          frb = InverseFourier[
              Abs[Fourier[dane[[fi]]]]
                  \^1.4]; 
			AppendTo[Wy, frb]]; 
		Wy]; \)\)

Funkcja detekcji maksimów funkcji autokorelacji i konwersja do Hz

\(\(\(DetMaksK = 
      Function[dane, 
		
        frange = {Floor[CzP\/CzMaks], Ceiling[CzP\/CzMin]}; 
		Wy = {}; 

        		For[fi = 1, fi <= Length[dane], \(fi++\), 
			fmaks = {0, 0}; 
          
			For[fj = frange[[1]], 
            fj <= frange[[2]], \(fj++\), 
            
				
            If[dane[[fi, fj]] > 
                  fmaks[[1]] && 
                dane[[fi, fj]] > 
                  dane[[fi, fj - 1]] && 
                dane[[fi, fj]] > 
                  dane[[fi, fj + 1]], 
              fmaks = {dane[[fi, fj]], 
                  fj}]]; 
			
          fpoly = \(-
              Fit[{dane[[fi, 
                    fmaks[[2]] - 1
                      ]], 
                  dane[[fi, 
                    fmaks[[2
                      ]]]], 
                  dane[[fi, 
                    fmaks[[2]] + 1
                      ]]}, {1, x, x\^2}, x]\); 
			
          fmaks[[2]] += 
            \((x /. \(FindMinimum[fpoly, {x, 1}]\)[[2
                    ]])\) - 2; 
			
          AppendTo[Wy, CzP/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

\(Timing[
	i = 0; 
	α = N[2\/π]; 
	Autok = {}; 

      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]; 
		
          If[wZP > 0, ZeroP = ZP[Okno, wZP], ZeroP = Okno]; 
		
          Widmo = OblWidma[ZeroP]; 
		
          WidPas = Podział[Widmo, Pasma]; 
		Xb = RekPasm[WidPas]; 
          
		Prost = Prostownik[Xb]; 
		Cb = FOAM[Prost, Pasma]; 
		
          Zb = N[α\ Xb + Cb]; 
		ZZb = ZeroPHamm[Zb]; 
		
          Rb = ACF[Zb]; 
		
          AppendTo[Autok, 
            ∑\+\(t = 1\)\%\(Length[Rb]\)Rb[[t
                ]]]]; 

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

      Autok = Abs[Autok]; ]\)

Detekcja maksimów i konwersja na Hz

FS = DetMaksK[Abs[Autok]];

Konwersja czasu

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

Wykres wyniku detekcji

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

Show[a, Graph];

Obliczenie błędu detekcji

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

Wykres błędu

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

Wartoć rednia miary błędu

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