Pitch tracking based on short-time processing

Include libraries

<< Miscellaneous`Audio`<< SignalProcessing`Master`

Signal Processing Function Definitions

Windowing

Windowing=Function[z,
		window=N[Table[DiscreteWindow[Blackman,WindowSize,winn],{winn, 0, WindowSize - 1}]];
		output = {};
		\(For[winx = 1, \ winx <= Length[window], \(winx++\), 
			
                output = 
                  Append[output, 
                    z[[winx]]*
                      window[[winx]]]]
                \);
		output];

Spectrum calculation

Spectrum = Function[specx, Abs[Fourier[Windowing[specx]]]];

Cepstrum calculation

\(\(Cepstrum = 
      Function[cepsx, 
		cepsy = Log[Spectrum[cepsx]]; 
		
        cepsy = Flatten[{
              Table[cepsy[[cepsn]], {
                  cepsn, \ WindowSize/2 + 1, \ WindowSize}], 
				
              Table[cepsy[[cepsn]], {
                  cepsn, \ WindowSize/2}]}]; 
		Abs[Fourier[cepsy]]
	]; 
    \)\)

Pitch tracking procedure

\(\(Track = 
      Function[trackdata, 
trackminceps = Floor[SmpRate/MaxPitch]; 

        trackmaxceps = Ceiling[SmpRate/MinPitch]; 
trackCepsPack = {}; 
		
        For[trackx = 1, trackx <= Length[trackdata] - WindowSize*Overlap, 
          trackx = trackx + WindowSize - WindowSize*Overlap, 
			
          trackblock = {}; 
			
          For[tracky = 0, tracky < WindowSize, \(tracky++\), 
				
            If[trackx + tracky > Length[trackdata], 
					
              trackblock = Append[trackblock, 0], 
					
              trackblock = 
                Append[trackblock, 
                  trackdata
                    [[trackx + tracky
                      ]]]]]; 
			
          trackCeps = Cepstrum[trackblock]; 
			
          trackCeps = 
            Table[trackCeps[[tracki]], {
                tracki, trackminceps, \ trackmaxceps}]; 
			
          trackCepsPack = Append[trackCepsPack, trackCeps]; 
		]; 
		\ 
        trackCepsPack
]; \)\)

Signal acquisition

SetDirectory["D:\\Mathematica\\Data\\"];

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

Additional Definitions

SmpRate=22050;
WindowSize=1024;
Overlap=(1/2);
MinPitch=200;
MaxPitch=2000;

Computations

Temporary calculations

res = Track[snd];

ListPlot3D[res,  PlotRange -> {5, 25}];