dawn@0: function [prd, fileinfo] = yin(x,p)
dawn@0: % YIN - Fundamental frequency (F0) of file or vector
dawn@0: %
dawn@0: % YIN('FILE') estimates and plots the F0 of signal in FILE.	 
dawn@0: %   F0 is plotted in octaves re: 440 Hz, together with residual 
dawn@0: %   (aperiodicity) and power measures.  A 'best' estimate is printed.
dawn@0: % YIN(X,SR) estimates and plots the F0 of array X assuming sampling rate SR (Hz).
dawn@0: %
dawn@0: % R=YIN(X) produces no plot but returns result in structure R:
dawn@0: %   R.f0: fundamental frequency in octaves re: 440 Hz
dawn@0: %   R.ap: aperiodicity measure (ratio of aperiodic to total power)
dawn@0: %   R.pwr: period-smoothed instantaneous power
dawn@0: %   (R also lists the parameters used by YIN)
dawn@0: %
dawn@0: % YIN(NAME,P) uses parameters stored in P:
dawn@0: %   P.minf0:    Hz - minimum expected F0 (default: 30 Hz)
dawn@0: %   P.maxf0:    Hz - maximum expected F0 (default: SR/(4*dsratio))
dawn@0: %   P.thresh:   threshold (default: 0.1)
dawn@0: %   P.relfag:   if ~0, thresh is relative to min of difference function (default: 1)
dawn@0: %   P.hop:      samples - interval between estimates (default: 32)
dawn@0: %   P.range:    samples - range of samples ([start stop]) to process
dawn@0: %   P.bufsize:  samples - size of computation buffer (default: 10000)
dawn@0: %	P.sr:		Hz - sampling rate (usually taken from file header)
dawn@0: %	P.wsize:	samples - integration window size (defaut: SR/minf0)
dawn@0: %	P.lpf:		Hz - intial low-pass filtering (default: SR/4)
dawn@0: %	P.shift		0: shift symmetric, 1: shift right, -1: shift left (default: 0)
dawn@0: %
dawn@0: % See 'yin.html' for more info.
dawn@0: % Version 28 July 2003.
dawn@0: 
dawn@0: % Alain de Cheveign? CNRS/Ircam, 2002.
dawn@0: % Copyright (c) 2002 Centre National de la Recherche Scientifique.
dawn@0: %
dawn@0: % Permission to use, copy, modify, and distribute this software without 
dawn@0: % fee is hereby granted FOR RESEARCH PURPOSES only, provided that this
dawn@0: % copyright notice appears in all copies and in all supporting 
dawn@0: % documentation, and that the software is not redistributed for any 
dawn@0: % fee (except for a nominal shipping charge). 
dawn@0: %
dawn@0: % For any other uses of this software, in original or modified form, 
dawn@0: % including but not limited to consulting, production or distribution
dawn@0: % in whole or in part, specific prior permission must be obtained from CNRS.
dawn@0: % Algorithms implemented by this software may be claimed by patents owned 
dawn@0: % by CNRS, France Telecom, Ircam or others.
dawn@0: %
dawn@0: % The CNRS makes no representations about the suitability of this 
dawn@0: % software for any purpose.  It is provided "as is" without express
dawn@0: % or implied warranty.
dawn@0: 
dawn@0: 
dawn@0: % Hidden parameters are integration window size (set equal to sr/minf0), search
dawn@0: % range for 'best neighboring estimate' (set equal to +/- sr/(2*minf0)), maximum
dawn@0: % expected width of period dip (stop/start ratio == 1.85), margin for "beam search"
dawn@0: % of final estimate (+1.8/-0.6 times the initial estimate).
dawn@0: 
dawn@0: 
dawn@0: % default parameter values ([]: to be determined)
dawn@0: minf0 = 30;			% Hz - minimum frequency
dawn@0: maxf0 = [];			% Hz - maximum frequency
dawn@0: wsize = []; 		% s - integration window size
dawn@0: lpf = [];			% Hz - lowpass prefiltering cutoff
dawn@0: thresh = 0.1;		% difference function threshold
dawn@0: relflag = 1;		% if true threshold is relative to global min of difference function
dawn@0: bufsize=10000;		% computation buffer size
dawn@0: hop = 32;			% samples - interval between estimates
dawn@0: range=[];			% range of file samples to process
dawn@0: sr=[];				% sampling rate
dawn@0: shift=0;			% flag to control the temporal shift of analysis windows (left/sym/right)
dawn@0: plotthreshold=0.2;	% aperiodicity above which plot is green or yellow
dawn@0: 
dawn@0: % if 2~=exist('allread')
dawn@0: % 	error('sf routines missing: put them in your path & try again');
dawn@0: % end
dawn@0: 
dawn@0: 
dawn@0: % handle parameters
dawn@0: if nargin<1; help yin; return; end
dawn@0: if nargin<2; p=[]; end
dawn@0: fileinfo=sf_info(x); if ~isempty(fileinfo.sr) p.sr=fileinfo.sr;	end % get sr from file
dawn@0: if fileinfo.nchans > 1
dawn@0: 	disp(['warning: using column 1 of ', num2str(fileinfo.nchans), '-column data']); 
dawn@0: end
dawn@0: if isa(p, 'double') p.sr=p; end
dawn@0: if ~isfield(p, 'sr'); p.sr=sr; end
dawn@0: if isempty(p.sr); error('YIN2: must specify SR'); end
dawn@0: if ~isfield(p, 'range') | isempty(p.range); p.range=[1 fileinfo.nsamples]; end
dawn@0: if ~isfield(p, 'minf0'); p.minf0=minf0; end
dawn@0: if ~isfield(p, 'thresh'); p.thresh=thresh; end
dawn@0: if ~isfield(p, 'relflag'); p.relflag=relflag; end
dawn@0: if ~isfield(p, 'bufsize'); p.bufsize=bufsize; end
dawn@0: if ~isfield(p, 'hop'); p.hop=hop; end
dawn@0: if ~isfield(p, 'maxf0'); p.maxf0=floor(p.sr/4); end % default
dawn@0: if ~isfield(p, 'wsize'); p.wsize=ceil(p.sr/p.minf0); end % default
dawn@0: if ~isfield(p, 'lpf'); p.lpf=p.sr/4; end % default
dawn@0: if mod(p.hop,1); error('hop should be integer'); end
dawn@0: if ~isfield(p, 'shift'); p.shift=shift; end % default
dawn@0: if ~isfield(p, 'plotthreshold'); p.plotthreshold=plotthreshold; end % default
dawn@0: 
dawn@0: % estimate period
dawn@0: r=yink(p,fileinfo);
dawn@0: prd=r.r1; % period in samples
dawn@0: ap0=r.r2; % gross aperiodicity measure
dawn@0: ap= r.r3; % fine aperiodicity measure
dawn@0: pwr=r.r4; % period-smoothed instantaneous power
dawn@0: f0 = log2(p.sr ./ prd) - log2(440); 	% convert to octaves re: 440 Hz
dawn@0: %figure();
dawn@0: %plot(prd);
dawn@0: %dlmwrite('pitch.txt',prd,'delimiter','\n');
dawn@0: %display(length(prd));
dawn@0: 
dawn@0: % load estimates and major parameters in result structure
dawn@0: clear r;
dawn@0: r.f0 = f0;
dawn@0: r.ap0 = ap0;
dawn@0: r.ap = ap;
dawn@0: r.pwr = pwr;
dawn@0: r.sr = p.sr;
dawn@0: r.range=p.range;
dawn@0: r.minf0 = p.minf0;
dawn@0: r.maxf0 = p.maxf0;
dawn@0: r.thresh=p.thresh;
dawn@0: r.relflag=p.relflag;
dawn@0: r.hop = p.hop;
dawn@0: r.bufsize = p.bufsize;
dawn@0: r.wsize = p.wsize;
dawn@0: r.lpf = p.lpf;
dawn@0: r.shift = p.shift;
dawn@0: r.plotthreshold=p.plotthreshold;
dawn@0: 
dawn@0: 
dawn@0: % plot estimates (if nargout == 0)
dawn@0: if nargout<1
dawn@0: 	if isnan(f0)
dawn@0: 		display('no estimates: signal too short or buffersize too small');
dawn@0: 		return;
dawn@0: 	end
dawn@0: 	% choose sample to report as "the" f0 of the entire signal
dawn@0: 	[mn, idx] = min(ap0);
dawn@0: 	best=f0(idx);
dawn@0: 	disp(['best: ', num2str(2^best*440), 'Hz (', note(best),...
dawn@0: 		') at ', num2str(idx/(p.sr/p.hop)), 's (aperiodic/total power: ', num2str(mn), ')']);
dawn@0: 	% plot f0 in 3 colors according to periodicity
dawn@0: 	good = f0;
dawn@0: 	good(find(ap0>p.plotthreshold*2)) = nan;
dawn@0: 	best = f0;
dawn@0: 	best(find(ap0>p.plotthreshold)) = nan;
dawn@0: 	subplot(211);
dawn@0: 	fsr=p.sr/p.hop;
dawn@0: 	nframes=size(prd,2);
dawn@0: 	if nframes <2; error('F0 track is too short to plot'); end
dawn@0: 	plot((1:nframes)/fsr, f0, 'y', (1:nframes)/fsr, good, 'g', (1:nframes)/fsr, best, 'b');
dawn@0: 	lo = max(min(f0),min(good)); hi=min(max(f0),max(good));
dawn@0: 	set(gca, 'ylim', [lo-0.5; hi+0.5]); 
dawn@0: 	set(gca, 'xlim', [1,nframes]/fsr);
dawn@0: 	set(get(gca,'ylabel'), 'string', 'Oct. re: 440 Hz');
dawn@0: 	set(get(gca,'xlabel'), 'string', 's');
dawn@0: 	% plot periodicity
dawn@0: 	ap0=max(0,ap0); ap=max(0,ap);
dawn@0: 	ap1=ap;
dawn@0: 	ap(find(ap0>plotthreshold)) = nan;
dawn@0: 	subplot(413); plot((1:nframes), ap0.^0.5, 'b'); 
dawn@0: 	%subplot(413); plot((1:nframes), ap0.^0.5, 'c', (1:nframes), ap1.^0.5, 'y',  (1:nframes), ap.^0.5, 'b'); 
dawn@0: 	set(gca, 'ylim', [0 1]); 
dawn@0: 	set(gca, 'xlim', [1, nframes]);
dawn@0: 	set(get(gca,'ylabel'), 'string', 'sqrt(apty)');
dawn@0: 	% plot power
dawn@0: 	subplot(414); plot((1:nframes), sqrt(pwr), 'b');	
dawn@0: 	set(gca, 'xlim', [1, nframes]);
dawn@0: 	set(get(gca,'ylabel'), 'string', 'sqrt(power)');
dawn@0: 	set(get(gca,'xlabel'), 'string', 'frames');
dawn@0:     if isa(x, 'double')
dawn@0:         set(gcf, 'Name', 'workspace matrix');
dawn@0:     else
dawn@0:         set (gcf, 'Name', x);
dawn@0:     end
dawn@0: end
dawn@0: 
dawn@0: 
dawn@0: % convert octave re 440 to note:
dawn@0: function s=note(o)
dawn@0: n=round(12*o);
dawn@0: cents = 100*(12*o-n);
dawn@0: oct=floor((n-3)/12)+5;
dawn@0: chroma=mod(n,12);
dawn@0: chromalist = {'A'; 'A#'; 'B'; 'C'; 'C#'; 'D'; 'D#'; 'E'; 'F'; 'F#';...
dawn@0: 	'G'; 'G#'};
dawn@0: cents = sprintf('%+.0f',cents);
dawn@0: s=[char(chromalist(chroma+1)),num2str(oct),' ',num2str(cents), ' cents'];