Re: Autocorrelation pitch detection ("Alain de Cheveigne'" )


Subject: Re: Autocorrelation pitch detection
From:    "Alain de Cheveigne'"  <alain.de.cheveigne@xxxxxxxx>
Date:    Fri, 3 Aug 2012 10:08:08 +0200
List-Archive:<http://lists.mcgill.ca/scripts/wa.exe?LIST=AUDITORY>

--Apple-Mail=_346A5641-60C7-4A70-9E23-6D49628ACCC4 Content-Transfer-Encoding: quoted-printable Content-Type: text/plain; charset=iso-8859-1 Hi Stuart, For the second you can use simple parabolic interpolation based on the 3 = samples nearest the peak. It's practically as accurate as since. Look = at my simple implementation of yin below. Interpolation is applied to = the difference function, but it works as well for the autocorrelation. For the first it is worth using a definition of autocorrelation that = does not introduce bias (rather than correcting the bias of the standard = definition). See discussion in my 1993 YIN paper. Alain On 2 Aug 2012, at 19:33, Stuart Rosen wrote: > Paul Boersma points out in his article 'ACCURATE SHORT-TERM ANALYSIS = OF THE FUNDAMENTAL FREQUENCY AND THE HARMONICS-TO-NOISE RATIO OF A = SAMPLED SOUND' = (http://www.fon.hum.uva.nl/paul/papers/Proceedings_1993.pdf) that there = are two ways to improve the use of an autocorrelation method to track = fundamental frequency in speech. One (correcting the autocorrelation for = the effects of the window function) is simple, but the interpolations he = suggests for best estimation of the lag and the height of that peak = (especially the sin(x)/x one) is not! Has anyone implemented these = procedures in Matlab by any chance, and are willing to share their code? >=20 > Thanks! >=20 > Yours - Stuart >=20 > --=20 > Stuart Rosen, PhD > Professor of Speech and Hearing Science > UCL Speech, Hearing and Phonetic Sciences > 2 Wakefield Street > London WC1N 1PF, England >=20 > Telephone numbers: > Office: (+ 44 [0]20) 7679 4077 > Admin: (+ 44 [0]20) 7679 4050 > Internal: ext 24077 > Fax: (+ 44 [0]20) 7679 4238 > Email: stuart@xxxxxxxx > Home page: http://www.phon.ucl.ac.uk/home/stuart >=20 > Request an inspection copy of Signals and Systems for Speech and = Hearing, 2nd edition by Stuart Rosen & Peter Howell > http://info.emeraldinsight.com/promo/signals.htm function r=3Dyin2(x,p); % YIN2: a simple implementation of the yin period-estimation algorithm % % yin2(x) : plot the period, power, and aperiodicity as a function of = time % % r=3Dyin2(x,p): use parameters in p, return result in r: % % r.prd: period=20 % r.ap: aperiodicity measure % =20 % p.maxprd: samples, maximum of search range [default 100] % p.minprd: samples, minimum of search range [default 2] % p.wsize: samples, window size [default maxprd] % p.hop: samples, frame period [default wsize] % p.thresh: threshold for period minimum [default: 0.1] % p.smooth: samples, size of low-pass smoothing window [default: = minprd/2] =20 if nargin>2; error('!'); end if nargin<2; p=3D[]; end if nargin<1; error('!'); end =20 % defaults if ~isfield(p, 'maxprd'); p.maxprd=3D100; end if ~isfield(p, 'minprd'); p.minprd=3D2; end if ~isfield(p, 'wsize'); p.wsize=3Dp.maxprd; end if ~isfield(p, 'hop'); p.hop=3Dp.wsize; end if ~isfield(p, 'thresh'); p.thresh=3D0.1; end if ~isfield(p, 'smooth'); p.smooth=3Dp.minprd/2; end =20 if min(size(x)) ~=3D 1; error('data should be 1D'); end x=3Dx(:); nsamples=3Dnumel(x); =20 nframes=3Dfloor((nsamples-p.maxprd-p.wsize)/p.hop); pwr=3Dzeros(1,nframes); prd=3Dzeros(1,nframes); ap=3Dzeros(1,nframes); =20 % shifted data x=3Dconvmtx(x,p.maxprd+1); x=3Dx(p.maxprd:end-p.maxprd,:); =20 for k=3D1:nframes =20 start=3D(k-1)*p.hop; % offset of frame xx=3Dx(start+1:start+p.wsize,:); d=3Dmean( (xx - repmat(xx(:,1),1,p.maxprd+1)).^2 )/2; % squared = difference function dd=3D d(2:end) ./ (cumsum(d(2:end)) ./ (1:(p.maxprd))); % = cumulative mean - normalized =20 % parabolic interpolation of all triplets to refine local minima min_pos=3D1:numel(dd); % nominal position of each sample x1=3Ddd(1:end-2); x2=3Ddd(2:end-1); x3=3Ddd(3:end); a=3D(x1+x3-2*x2)/2; b=3D(x3-x1)/2; shift=3D-b./(2*a); % offset of interpolated minimum re = current sample val=3Dx2-b.^2./(4*a); % value of interpolated minimum =20 % replace all local minima by their interpolated value,=20 idx=3D 1 + find(x2<x1 & x2<x3); dd(idx)=3Dval(idx-1); min_pos(idx)=3Dmin_pos(idx-1)+shift(idx-1); =20 % find index of first min below threshold a=3Ddd<p.thresh; if isempty(find(a)) [~,prd0]=3Dmin(dd); % none below threshold, take global min = instead else b=3Dmin(find(a)); % left edge c=3Dmin(b*2,numel(a)); [~,prd0]=3Dmin(dd(b:(c-1)));=20 prd0=3Db+prd0-1; end =20 prd=3Dmin_pos(prd0)+1; =20 if prd>2 & prd<numel(dd) & d(prd0)<d(prd0-1) & d(prd0)<d(prd0+1) =20 % refine by parabolic interpolation of raw difference function=20= x1=3Dd(prd-1); x2=3Dd(prd); x3=3Dd(prd+1); a=3D(x1+x3-2*x2)/2; b=3D(x3-x1)/2; shift=3D-b./(2*a); % offset of interpolated minimum re = current sample val=3Dx2-b.^2./(4*a); % value of interpolated minimum prd=3Dprd+shift-1; end =20 % aperiodicity frac=3Dprd-floor(prd); if frac=3D=3D0 yy=3Dxx(:,prd); else yy=3D(1-frac)*xx(:,floor(prd+1))+frac*xx(:,floor(prd+1)+1); % = linear interpolation end pwr=3D(mean(xx(:,1).^2) + mean(yy.^2))/2; % average power over fixed = and shifted windows res=3Dmean(((xx(:,1) - yy)).^2) / 2; ap=3Dres/pwr; =20 =20 r.prd(k)=3Dprd; r.ap(k)=3Dap; =20 end =20 if nargout=3D=3D0;=20 subplot 211; plot(r.prd); title('period'); xlabel('frame'); = ylabel('samples'); subplot 212; plot(r.ap); title('periodicity'); xlabel('frame'); r=3D[];=20 end =20 --Apple-Mail=_346A5641-60C7-4A70-9E23-6D49628ACCC4 Content-Transfer-Encoding: quoted-printable Content-Type: text/html; charset=iso-8859-1 <html><head></head><body style=3D"word-wrap: break-word; = -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi = Stuart,<div><br></div><div>For the second you can use simple parabolic = interpolation based on the 3 samples nearest the peak. &nbsp;It's = practically as accurate as since. &nbsp;Look at my simple implementation = of yin below. &nbsp;Interpolation is applied to the difference function, = but it works as well for the = autocorrelation.</div><div><br></div><div>For the first it is worth = using a definition of autocorrelation that does not introduce bias = (rather than correcting the bias of the standard definition). &nbsp;See = discussion in my 1993 YIN = paper.</div><div><br></div><div>Alain</div><div><br><div><div>On 2 Aug = 2012, at 19:33, Stuart Rosen wrote:</div><br = class=3D"Apple-interchange-newline"><blockquote type=3D"cite"> =20 <meta http-equiv=3D"content-type" content=3D"text/html; = charset=3DISO-8859-1"> =20 <div bgcolor=3D"#FFFFFF" text=3D"#000000"> Paul Boersma points out in his article 'ACCURATE SHORT-TERM ANALYSIS OF THE FUNDAMENTAL FREQUENCY AND THE HARMONICS-TO-NOISE RATIO OF A SAMPLED SOUND' (<a class=3D"moz-txt-link-freetext" = href=3D"http://www.fon.hum.uva.nl/paul/papers/Proceedings_1993.pdf">http:/= /www.fon.hum.uva.nl/paul/papers/Proceedings_1993.pdf</a>) that there are two ways to improve the use of an autocorrelation method to track fundamental frequency in speech. One (correcting the autocorrelation for the effects of the window function) is simple, but the interpolations he suggests for best estimation of the lag and the height of that peak (especially the sin(x)/x one) is = not!&nbsp; Has anyone implemented these procedures in Matlab by any chance, and are willing to share their code?<br> <br> Thanks!<br> <br> Yours - Stuart<br> <br> <div class=3D"moz-signature">-- <br> <hr> Stuart Rosen, PhD<br> Professor of Speech and Hearing Science<br> UCL Speech, Hearing and Phonetic Sciences<br> 2 Wakefield Street<br> London WC1N 1PF, England<br> <br> Telephone numbers: <ul> <li>Office: (+ 44 [0]20) 7679 4077</li> <li>Admin: (+ 44 [0]20) 7679 4050</li> <li>Internal: ext 24077</li> <li>Fax: (+ 44 [0]20) 7679 4238</li> </ul> Email: <a class=3D"moz-txt-link-abbreviated" = href=3D"mailto:stuart@xxxxxxxx">stuart@xxxxxxxx</a><br> Home page: <a class=3D"moz-txt-link-freetext" = href=3D"http://www.phon.ucl.ac.uk/home/stuart">http://www.phon.ucl.ac.uk/h= ome/stuart</a><br> <br> Request an inspection copy of <i>Signals and Systems for Speech and Hearing, 2nd edition</i> by Stuart Rosen &amp; Peter Howell<br> <a class=3D"moz-txt-link-freetext" = href=3D"http://info.emeraldinsight.com/promo/signals.htm">http://info.emer= aldinsight.com/promo/signals.htm</a><br> <hr> </div> </div> </blockquote></div><br></div><div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">function</span> r=3Dyin2(x,p);</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; color: rgb(30, 152, 65); ">% = YIN2: a simple implementation of the yin period-estimation = algorithm</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(30, 152, 65); ">%&nbsp; yin2(x) : plot the period, power, and = aperiodicity as a function of time</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); = ">%</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; r=3Dyin2(x,p): = use parameters in p, return result in r:</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; color: rgb(30, 152, 65); = ">%</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; &nbsp; r.prd: = period&nbsp;</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; &nbsp; r.ap: = aperiodicity measure</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">% &nbsp;</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(30, 152, 65); ">%&nbsp; &nbsp; p.maxprd: samples, maximum of search = range [default 100]</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; &nbsp; p.minprd: = samples, minimum of search range [default 2]</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(30, 152, 65); ">%&nbsp; &nbsp; p.wsize: samples, window size = [default maxprd]</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; &nbsp; p.hop: = samples, frame period [default wsize]</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; color: rgb(30, 152, 65); = ">%&nbsp; &nbsp; p.thresh: threshold for period minimum [default: = 0.1]</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); ">%&nbsp; &nbsp; p.smooth: = samples, size of low-pass smoothing window [default: minprd/2]</div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: = #1e9841; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">if</span> nargin&gt;2; error(<span style=3D"color: = #b34de9">'!'</span>); <span style=3D"color: = #063ff4">end</span></div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; "><span style=3D"color: #063ff4">if</span> = nargin&lt;2; p=3D[]; <span style=3D"color: #063ff4">end</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = "><span style=3D"color: #063ff4">if</span> nargin&lt;1; error(<span = style=3D"color: #b34de9">'!'</span>); <span style=3D"color: = #063ff4">end</span></div><p style=3D"margin: 0.0px 0.0px 0.0px 0.0px; = font: 10.0px Courier; color: #063ff4; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); ">% = defaults</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; "><span style=3D"color: #063ff4">if</span> = ~isfield(p, <span style=3D"color: #b34de9">'maxprd'</span>); = p.maxprd=3D100; <span style=3D"color: #063ff4">end</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = "><span style=3D"color: #063ff4">if</span> ~isfield(p, <span = style=3D"color: #b34de9">'minprd'</span>); p.minprd=3D2; <span = style=3D"color: #063ff4">end</span></div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">if</span> ~isfield(p, <span style=3D"color: = #b34de9">'wsize'</span>); p.wsize=3Dp.maxprd; <span style=3D"color: = #063ff4">end</span></div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; "><span style=3D"color: #063ff4">if</span> = ~isfield(p, <span style=3D"color: #b34de9">'hop'</span>); p.hop=3Dp.wsize;= <span style=3D"color: #063ff4">end</span></div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">if</span> ~isfield(p, <span style=3D"color: = #b34de9">'thresh'</span>); p.thresh=3D0.1; <span style=3D"color: = #063ff4">end</span></div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; "><span style=3D"color: #063ff4">if</span> = ~isfield(p, <span style=3D"color: #b34de9">'smooth'</span>); = p.smooth=3Dp.minprd/2; <span style=3D"color: #063ff4">end</span></div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: = #063ff4; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">if</span> min(size(x)) ~=3D 1; error(<span style=3D"color: = #b34de9">'data should be 1D'</span>); <span style=3D"color: = #063ff4">end</span></div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">x=3Dx(:);</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">nsamples=3Dnumel(x);</div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; = min-height: 12.0px">&nbsp;<br class=3D"webkit-block-placeholder"></p><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">nframes=3Dfloor((nsamples-p.maxprd-p.wsize)/p.hop);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">pwr=3Dzeros(1,nframes);</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">prd=3Dzeros(1,nframes);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">ap=3Dzeros(1,nframes);</div><p style=3D"margin: 0.0px 0.0px 0.0px = 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); ">% shifted = data</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">x=3Dconvmtx(x,p.maxprd+1);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">x=3Dx(p.maxprd:end-p.maxprd,:);</div><p style=3D"margin: 0.0px 0.0px = 0.0px 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">for</span> k=3D1:nframes</div><p style=3D"margin: 0.0px 0.0px = 0.0px 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; start=3D(k-1)*p.hop; = <span style=3D"color: #1e9841">% offset of frame</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; xx=3Dx(start+1:start+p.wsize,:);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; d=3Dmean( (xx - repmat(xx(:,1),1,p.maxprd+1)).^2 )/2; = &nbsp; &nbsp; <span style=3D"color: #1e9841">% squared difference = function</span></div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; dd=3D d(2:end) ./ (cumsum(d(2:end)) = ./ (1:(p.maxprd))); &nbsp; <span style=3D"color: #1e9841">% cumulative = mean - normalized</span></div><p style=3D"margin: 0.0px 0.0px 0.0px = 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp;&nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; </span>% parabolic interpolation = of all triplets to refine local minima</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; color: rgb(30, 152, 65); = "><span style=3D"color: #000000">&nbsp; &nbsp; = min_pos=3D1:numel(dd);&nbsp; &nbsp; </span>% nominal position of each = sample</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; x1=3Ddd(1:end-2);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; x2=3Ddd(2:end-1);</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; = x3=3Ddd(3:end);</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; a=3D(x1+x3-2*x2)/2;</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; b=3D(x3-x1)/2;</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; shift=3D-b./(2*a);&nbsp; &nbsp; = &nbsp; &nbsp; </span>% offset of interpolated minimum re current = sample</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); "><span style=3D"color: = #000000">&nbsp; &nbsp; val=3Dx2-b.^2./(4*a); &nbsp; &nbsp; </span>% = value of interpolated minimum</div><p style=3D"margin: 0.0px 0.0px 0.0px = 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp;&nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; </span>% replace all local minima = by their interpolated value,&nbsp;</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; idx=3D 1 + = find(x2&lt;x1 &amp; x2&lt;x3);</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; = dd(idx)=3Dval(idx-1);</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; = min_pos(idx)=3Dmin_pos(idx-1)+shift(idx-1);</div><p style=3D"margin: = 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; min-height: = 12.0px">&nbsp;&nbsp; &nbsp;<br class=3D"webkit-block-placeholder"></p><div= style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(30, 152, 65); "><span style=3D"color: #000000">&nbsp; &nbsp; = </span>% find index of first min below threshold</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; a=3Ddd&lt;p.thresh;</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; <span style=3D"color: = #063ff4">if</span> isempty(find(a))</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; &nbsp; &nbsp; [~,prd0]=3Dmin(dd); = </span>% none below threshold, take global min instead</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; <span style=3D"color: #063ff4">else</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; &nbsp; &nbsp; b=3Dmin(find(a)); <span style=3D"color: = #1e9841">% left edge</span></div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = c=3Dmin(b*2,numel(a));</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = [~,prd0]=3Dmin(dd(b:(c-1)));&nbsp;</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = prd0=3Db+prd0-1;</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; <span style=3D"color: = #063ff4">end</span></div><p style=3D"margin: 0.0px 0.0px 0.0px 0.0px; = font: 10.0px Courier; min-height: 12.0px">&nbsp;&nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; = prd=3Dmin_pos(prd0)+1;</div><p style=3D"margin: 0.0px 0.0px 0.0px 0.0px; = font: 10.0px Courier; min-height: 12.0px">&nbsp;&nbsp; &nbsp; &nbsp; = &nbsp;<br class=3D"webkit-block-placeholder"></p><div style=3D"margin-top:= 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; ">&nbsp; &nbsp; <span = style=3D"color: #063ff4">if</span> prd&gt;2 &amp; prd&lt;numel(dd) &amp; = d(prd0)&lt;d(prd0-1) &amp; d(prd0)&lt;d(prd0+1)</div><p style=3D"margin: = 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; min-height: = 12.0px">&nbsp;&nbsp; &nbsp; &nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; &nbsp; &nbsp; </span>% refine by = parabolic interpolation of raw difference function&nbsp;</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; &nbsp; &nbsp; x1=3Dd(prd-1);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; &nbsp; &nbsp; x2=3Dd(prd);</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = x3=3Dd(prd+1);</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = a=3D(x1+x3-2*x2)/2;</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = b=3D(x3-x1)/2;</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(30, 152, 65); "><span style=3D"color: = #000000">&nbsp; &nbsp; &nbsp; &nbsp; shift=3D-b./(2*a);&nbsp; &nbsp; = &nbsp; &nbsp; </span>% offset of interpolated minimum re current = sample</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; val=3Dx2-b.^2./(4*a); = &nbsp; &nbsp; <span style=3D"color: #1e9841">% value of interpolated = minimum</span></div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; &nbsp; &nbsp; = prd=3Dprd+shift-1;</div><div style=3D"margin-top: 0px; margin-right: = 0px; margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; <span style=3D"color: = #063ff4">end</span></div><p style=3D"margin: 0.0px 0.0px 0.0px 0.0px; = font: 10.0px Courier; color: #063ff4; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; color: rgb(30, 152, 65); "><span = style=3D"color: #000000">&nbsp; &nbsp; </span>% aperiodicity</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; frac=3Dprd-floor(prd);</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; ">&nbsp; &nbsp; <span = style=3D"color: #063ff4">if</span> frac=3D=3D0</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; &nbsp; &nbsp; yy=3Dxx(:,prd);</div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; <span style=3D"color: #063ff4">else</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; &nbsp; &nbsp; = yy=3D(1-frac)*xx(:,floor(prd+1))+frac*xx(:,floor(prd+1)+1); <span = style=3D"color: #1e9841">% linear interpolation</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; <span style=3D"color: #063ff4">end</span></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(30, 152, 65); "><span style=3D"color: #000000">&nbsp; &nbsp; = pwr=3D(mean(xx(:,1).^2) + mean(yy.^2))/2; </span>% average power over = fixed and shifted windows</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; res=3Dmean(((xx(:,1) = - yy)).^2) / 2;</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; ">&nbsp; &nbsp; ap=3Dres/pwr;</div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; = min-height: 12.0px">&nbsp;&nbsp; &nbsp;<br = class=3D"webkit-block-placeholder"></p><p style=3D"margin: 0.0px 0.0px = 0.0px 0.0px; font: 10.0px Courier; min-height: 12.0px">&nbsp;&nbsp; = &nbsp; &nbsp; &nbsp;<br class=3D"webkit-block-placeholder"></p><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; = ">&nbsp; &nbsp; r.prd(k)=3Dprd;</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; r.ap(k)=3Dap;</div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; = min-height: 12.0px">&nbsp;<br class=3D"webkit-block-placeholder"></p><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 10px/normal Courier; color: = rgb(6, 63, 244); ">end</div><p style=3D"margin: 0.0px 0.0px 0.0px 0.0px; = font: 10.0px Courier; color: #063ff4; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; "><span style=3D"color: = #063ff4">if</span> nargout=3D=3D0;&nbsp;</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; ">&nbsp; &nbsp; subplot <span = style=3D"color: #b34de9">211</span>; plot(r.prd); title(<span = style=3D"color: #b34de9">'period'</span>); xlabel(<span style=3D"color: = #b34de9">'frame'</span>); ylabel(<span style=3D"color: = #b34de9">'samples'</span>);</div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 10px/normal Courier; ">&nbsp; &nbsp; subplot <span = style=3D"color: #b34de9">212</span>; plot(r.ap); title(<span = style=3D"color: #b34de9">'periodicity'</span>); xlabel(<span = style=3D"color: #b34de9">'frame'</span>);</div><div style=3D"margin-top: = 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: = normal normal normal 10px/normal Courier; ">&nbsp; &nbsp; = r=3D[];&nbsp;</div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 10px/normal Courier; color: rgb(6, 63, 244); ">end</div><p = style=3D"margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: = #063ff4; min-height: 12.0px">&nbsp;<br = class=3D"webkit-block-placeholder"></p><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 12px/normal Courier; min-height: 14px; = "><br></div></div><div style=3D"margin-top: 0px; margin-right: 0px; = margin-bottom: 0px; margin-left: 0px; font: normal normal normal = 12px/normal Courier; min-height: 14px; "><br></div><div = style=3D"margin-top: 0px; margin-right: 0px; margin-bottom: 0px; = margin-left: 0px; font: normal normal normal 12px/normal Courier; = min-height: 14px; "><br></div><div style=3D"margin-top: 0px; = margin-right: 0px; margin-bottom: 0px; margin-left: 0px; font: normal = normal normal 12px/normal Courier; min-height: 14px; = "><br></div></body></html>= --Apple-Mail=_346A5641-60C7-4A70-9E23-6D49628ACCC4--


This message came from the mail archive
/var/www/postings/2012/
maintained by:
DAn Ellis <dpwe@ee.columbia.edu>
Electrical Engineering Dept., Columbia University