OK, Robert, just for you... I am posting it a second time.
This is cut-and-paste Matlab routine - just fix the wrapped around
comment lines.
Add some signal pre-processing (like low-pass filtering) as you wish...
You are still confused, my friend, about pitch vs. fundamental
frequency (or, rather, fundamental period)
Pitch is an auditory concept, it's tied to human perception.
PDA is a mathematical method measuring fundamental period (for
correlation, AMDF, ASDF, or my method) or fundamental frequency (for
Fourier-based methods) of a periodic signal.
(After "true" f0 is measured an adjustment can be made to accomodate
for limitations of a human auditory perception)
It is that simple: perception vs. strict math
And, BTW, I am not trying to sell you anything: I am a scientist and
inventor, not a salesman, for your information
Regards
**************************************************************************
This is a prototype software implementation of the invention.
It illustrates basic principles and operation of the new methods.
The m-files below (ptrack.m, pframe.m, embed.m, xdist.m, nhistmax.m,
histsmooth.m, localpeaks.m and, optionally, svdembed.m) should be
copied
to Matlab work directory.
Copyright (C) 2001 Dmitry E. Terez
DISCLAIMER: This software is provided WITHOUT ANY WARRANTY, without
even
implied statement of merchantability or fitness for a particular
purpose.
******************************** ptrack.m
********************************
function ptr =3D ptrack(x)
%PTRACK Compute pitch period track.
% This MATLAB function computes pitch period track and performs
% voiced/unvoiced segmentation of a speech signal.
%
% x - input speech signal vector (MATLAB double-precision array).
% ptr - output vector of pitch period estimates for all fixed-size
frames.
% Each frame is assigned either a pitch period value in integer
samples
% (for frames determined to be voiced) or zero (for unvoiced
frames).
% To obtain fundamental frequency f0 values from pitch period
values,
% one needs to compute Fs/ptr for each of the voiced frames, where
% Fs is a sampling frequency in Hz and ptr is pitch period in
samples.
% (For unvoiced frames f0 is assigned zero value).
% Further smoothing of the raw output pitch tracks may be desirable.
%
% This program was tested under MATLAB 5.3.1 on Windows PC platform.
% No additional toolboxes are required.
% Parameters are approximately chosen to produce correct results in
% most cases on clean female speech and male speech with f0 > 100Hz
% from TIMIT database (16 kHz sampling, 16 bits per sample).
% For other sampling rates and speech types (e.g. male speech with
f0<100Hz)
% some parameters (e.g. frame size) will need to be modified.
% Frame parameters:
nwin =3D 200; % frame size in samples
nstep =3D 100; % frame step in samples
% Pitch search range is 100 Hz - 400 Hz (at 16 kHz sampling rate):
low =3D 40; % smallest pitch period to search for (in samples)
high =3D 160; % largest pitch period to search for (in samples)
% Parameters for time-delay embedding:
ndim =3D 3; % embedding dimension
ndel =3D 10; % time delay in samples
% Parameters for pitch tracking and voicing decision:
% emin is a minimum RMS energy of voiced frame.
% This value may require adjustment.
% The present value is approximately chosen for clean speech
% signals with an amplitude normalized to be in [-1 +1] range
emin =3D 0.015;
pdiff =3D 6; % max. allowed pitch period difference in samples
% between adjacent voiced frames
pdiff2 =3D 2.0*pdiff;
nback =3D 5; % max. number of frames to track backward
ne =3D nwin - (ndim-1)*ndel; % Number of vectors in state space
% after embedding a speech frame
% Global variable KK is an array with 2 columns storing all
% possible non-repeating combinations of 2 integers
% with their absolute difference between low and high.
% It is computed only once using Matlab 5 function nchoosek
% and subsequently reused for each frame computation.
np =3D 1:ne;
nn =3D nchoosek(np,2);
dt =3D nn(:,2)-nn(:,1);
isel =3D find(dt > low & dt < high);
global KK
KK =3D nn(isel,:);
% Determine input length and required number of frames
nx=3Dlength(x);
nframes=3Dfloor(1+(nx-nwin)/nstep);
x=3Dx(:); % Make input signal vector a column
% Normalize input signal amplitude to be in [-1 1] range
maxval=3Dmax(abs(x));
x=3Dx/maxval;
% Initialize arrays
ptr =3D zeros(1,nframes); % Output vector to store final pitch values
pcan =3D cell(1,nframes); % Cell array to store pitch candidates for
all frames
rel =3D zeros(1,nframes); % Vector to store reliability measures for
all frames
en =3D zeros(1,nframes); % Vector to store RMS energy for all frames
for i=3D1:nframes % The main loop going through all frames
nstart=3D(i-1)*nstep+1;
nend=3Dnstart+nwin-1;
xf =3D x(nstart:nend); % current effective frame
% Process current effective frame
[pe, rm, nsel, rsel, xen] =3D pframe(xf, ndim, ndel, low, high);
pcan{i} =3D pe; % array of pitch candidates
rel(i) =3D rm; % reliability measure
en(i) =3D xen; % RMS frame energy
% Track pitch delaying final decisions by one or more frames.
if(i > 2) % At least previous frame and next frame are needed
jp =3D i-2; % Previous frame
jc =3D i-1; % Current frame (delayed by one frame)
jn =3D i; % Next frame (same as current effective frame)
if (en(jc) < emin) % If energy is too low assign unvoiced value
ptr(jc) =3D 0;
else % Energy is sufficient
if (rel(jc) =3D=3D 1) % current frame is reliable
if (ptr(jp) > 0) % previous frame was found voiced
if(abs(pcan{jc}(1)-ptr(jp))<pdiff) %cur. consistent with
prev.
ptr(jc) =3D pcan{jc}(1); % assign pitch (lowest multiple)
% Otherwise check next frame to see if it's reliable and
% consistent with previous frame (current frame fall-out)
elseif(rel(jn)=3D=3D1 & abs(pcan{jn}(1)-ptr(jp)) < pdiff2)
ptr(jc)=3Dround(0.5*(ptr(jp)+pcan{jn}(1)));
% Check if next frame is reliable and consistent with current
elseif (rel(jn)=3D=3D1 & abs(pcan{jn}(1)-pcan{jc}(1)) < pdiff)
ptr(jc)=3Dpcan{jc}(1); %can be real pitch doubling or
halfing
else % If all checks fail assign unvoiced value to current
ptr(jc) =3D 0;
end
else % previous frame was found unvoiced
% Check for a start of voiced segment: need two consecutive,
% reliable and consistent with each other frames
if(rel(jn) =3D=3D 1 & abs(pcan{jn}(1)-pcan{jc}(1)) < pdiff)
ptr(jc) =3D pcan{jc}(1);
% If start of voiced segment is determined, try to track
% pitch period backward for up to nback frames
jj =3D jc-1;
found =3D 1;
nf =3D 0;
while ( jj > 0 & ptr(jj) =3D=3D 0 & en(jj) > emin ...
& found =3D=3D 1 & nf < nback)
parr =3D pcan{jj}; % vector of pitch candidates
[mindiff,k] =3D min(abs(parr - ptr(jj+1)));
if (mindiff < pdiff) % If found a match
ptr(jj) =3D parr(k); % then assign the value found
else % Stop backward tracking
found =3D 0;
end
jj =3D jj - 1;
nf =3D nf + 1;
end
else % Otherwise assign unvoiced value for now
ptr(jc) =3D 0;
end
end
else % current frame is unreliable
if (ptr(jp)>0) % previous frame was found voiced
% Search current frame for best match to previous
parr =3D pcan{jc}; % vector of current pitch candidates
[mindiff, k] =3D min(abs(parr - ptr(jp)));
if (mindiff < pdiff) % If found a good match
ptr(jc) =3D parr(k); % then assign the value found
else % Otherwise check if next frame
if(en(jn)>emin & rel(jn)=3D=3D1 & ... % is reliable
abs(pcan{jn}(1)-ptr(jp))<pdiff2) % and matches prev.
ptr(jc)=3Dround(0.5*(ptr(jp)+pcan{jn}(1)));
else % If all checks fail assign unvoiced value
ptr(jc) =3D 0;
end
end
else % previous frame was found unvoiced
ptr(jc) =3D 0; % assign unvoiced value to current frame
end
end % current frame reliable or unreliable
end % if energy < emin
end % if i > 2
end % Main for loop
**************************************************************************
***************************** pframe.m
***********************************
function [pcan,rel,nsel,rsel,xen] =3D pframe(xf,ndim,ndel,low,high)
%PFRAME Determine pitch candidates for a speech frame.
% This function determines pitch period candidates for
% a signal frame of some predetermined number of samples.
%
% INPUT:
% xf - input signal vector (speech frame)
% ndim - embedding dimension
% ndel - delay value for time-delay embedding (in integer samples)
% low - lower bound for pitch search in integer samples
% high - upper bound for pitch search in integer samples
%
% OUTPUT:
% pcan - vector of pitch period candidates (in integer samples),
% ordered by value from lowest to highest (if rel=3D=3D1),
% or by corresponding peak magnitude from highest peak
% to lowest (if rel=3D=3D0).
% rel - frame reliability measure:
% rel=3D1 if the number of pitch candidates <=3D 3
% and they are integer multiples (with some tolerance),
% or there is only one pitch candidate;
% rel=3D0 otherwise (or if no candidates are found).
% nsel - number of distances selected for a histogram
% rsel - neigborhood radius
% xen - RMS frame energy
% Parameters (in relative units - frame size independent):
rmax =3D 0.15; % maximal neigborhood radius in state space
hmax =3D 0.9; % high reference peak magnitude
href =3D 0.8; % reference peak magnitude
hmin =3D 0.7; % low reference peak magnitude
hrelmin =3D 0.6; % minimal magnitude of highest peak to have rel=3D1
mtol =3D 0.05; % relative tolerance for testing pitch multiples
thld =3D 0.4; % threshold level
minsel =3D 2; % min. number of distances as fraction of nwin
maxpeaks =3D 6; % max. number of pitch candidates to retain
smooth =3D 1; % perform histogram smoothing before peak searching
minsep =3D 8; % min. separation allowed between peaks (in samples)
nx =3D length(xf); % number of samples in input frame
nselmin =3D minsel*nx; % min. number of distances to select for
histogram
% Embed signal using time-delay embedding
% (Normalization is also performed here in order for all
% output vectors to fit into a unit cube in state space)
xe =3D embed(xf, ndim, ndel);
% To use alternative SVD-embedding uncomment next 2 lines
%nsvd =3D 20; % window size for SVD-embedding (can be adjusted)
%xe =3D svdembed(xf, ndim, nsvd);
% nxe is the number of vectors in state space
[nxe, junk] =3D size(xe);
maxdt =3D nxe - 1; % maximal time separation (in samples) between
% vectors in state space (histogram length)
% Find and compute all spatial distances less than rmax
% and corresponding temporal separations
xd =3D xdist(xe, rmax);
% ntotal is the total number of distances less than rmax
[ntotal, junk] =3D size(xd);
dt =3D xd(:,2)'; % Extract vector of temporal separations
% (sorted by corresponding spatial distances
% from smallest to largest)
% First, select all found distances closer than rmax,
% compute normalized histogram and find highest peak m0
nsel =3D ntotal;
dtsel =3D dt(1:nsel);
[nhist, m0] =3D nhistmax(dtsel, maxdt, low, high);
% Try to adjust neigborhood radius in order to bring
% the highest peak to some desirable range, if possible
if ( m0 > hmax & ntotal > nselmin )
nsel =3D nselmin;
dtsel =3D dt(1:nsel);
[nhist, m0] =3D nhistmax(dtsel, maxdt, low, high);
if ( m0 < hmin )
nsel =3D min(floor(nsel*href/m0),ntotal);
dtsel =3D dt(1:nsel);
[nhist, m0] =3D nhistmax(dtsel, maxdt, low, high);
end
end
hist =3D nhist-thld*m0; % Threshold histogram
hist1=3Dhist.*(sign(hist)+1)./2; % and half-way rectify
% (Optionally) smooth histogram before searching for peaks
if (smooth =3D=3D 1)
hist1 =3D histsmooth(hist1);
end
% Find all local peaks between low and high bounds
% (with optional constraint on their proximity to each other)
peaks =3D localpeaks(hist1,low,high,minsep);
[npeaks, junk] =3D size(peaks); % total number of peaks found
% Analyze number and positions of found peaks
if (npeaks =3D=3D 0) % If no peaks are found
pcan =3D 0; % assign unvoiced value
rel =3D 0; % frame is unreliable
elseif (npeaks > 0 & npeaks <=3D 3) % 1, 2 or 3 peaks are found
pi =3D sort(peaks(:,2)); % Sort peaks by position
% Check for multiples (with some relative tolerance mtol)
if (npeaks =3D=3D 1)
mult =3D 1;
elseif (npeaks =3D=3D 2 & abs(pi(2)-2.0*pi(1))/pi(1) < mtol)
mult =3D 1;
elseif (npeaks =3D=3D 3 & abs(pi(2)-2.0*pi(1))/pi(1) < mtol &...
abs(pi(3)-3.0*pi(1))/pi(1) < mtol)
mult =3D 1;
else
mult =3D 0;
end
if (mult =3D=3D 1 & m0 > hrelmin) % if multiples and mo>hrelmin
rel =3D 1; % proclaim current frame to be reliable
else
rel =3D 0; % otherwise frame is unreliable
end
pcan =3D pi; % pitch period candidates
else % more than 3 peaks found
np =3D min(npeaks, maxpeaks);
pcan =3D peaks(1:np,2); % retain np highest peaks as pitch cand.
rel =3D 0; % frame is unreliable
end % npeaks
rsel =3D xd(nsel,1); % chosen neighborhood radius
xen =3D sqrt(sum(xf.*xf)/nx); % RMS frame energy
**************************************************************************
****************************** embed.m
***********************************
function xe =3D embed(x, ndim, ndel)
%EMBED Embed a signal using time-delay embedding.
% x - input signal vector
% ndim - embedding dimension
% ndel - time delay in samples
% xe - output matrix with ndim columns and nr rows
% (Each column represents a separate dimension and
% each row represents a vector in state space)
nx =3D length(x); % input length
% Normalize input vector x to be in [0 1] range in order
% for a reconstructed trajectory to fit into a unit cube
% (Alternatively, this can be done to each column of xe)
xmax =3D max(x);
xmin =3D min(x);
xrange =3D xmax-xmin;
if (xrange > 0)
x =3D (x - xmin)/xrange;
else
x =3D zeros(size(x));
end
% Calculate the number of rows in the output array xe
nr =3D nx - (ndim-1)*ndel; % Number of vectors in state space
xe=3Dzeros(nr,ndim); % Preallocate space for output array
% Embed using method of delays
for i =3D 1:ndim
xe(:,i) =3D x(1+ndel*(i-1) : nr+ndel*(i-1));
end
**************************************************************************
****************************** xdist.m
***********************************
function xd =3D xdist(xe, rmax)
%XDIST Select and compute distances in state space and time
separations.
% This function finds all distances less than rmax between pairs of
% vectors in state space and computes corresponding time
separations.
% INPUT:
% xe - input array with embedded signal (vectors in state space)
% rmax - chosen neighborhood radius in state space
% OUTPUT:
% xd - 2-column array of computed spatio-temporal distances:
% Col. 1 has Euclidean distances less than rmax;
% Col. 2 has corresponding time separations (in samples).
% Array xd is sorted in ascending order on the 1st column.
%
% (This function can be replaced with any suitable fast
% neighbor-searching procedure to speed-up computations.)
global KK % Global array KK of possible index combinations is reused
rmax2 =3D rmax*rmax; % Square of max. Euclidean distance
[mx,ndim] =3D size(xe); % mx - number of vectors in state space
% ndim - embedding dimension
[mk,junk] =3D size(KK); % mk - number of vector pairs
% Initialize a big array of spatio-temporal distances
% and fill 2nd column with temporal separations
xd =3D [zeros(mk,1) (KK(:,2)-KK(:,1))];
% Compute all squared Euclidean distances and fill 1st column
% (Distance measures other than Euclidean are also possible)
for i=3D1:ndim
xd(:,1)=3Dxd(:,1) + (xe(KK(:,2),i)-xe(KK(:,1),i)).^2;
end
% Select only squared Euclidean distances less than rmax2
isel =3D find( xd(:,1) < rmax2 );
xd =3D xd(isel,:);
% Sort in ascending order based on distances in the first column
xd =3D sortrows(xd);
% For computational efficiency squared values of Euclidean distances
% should be used instead of actual distances:
% square root operation is not needed and is retained here
% only for illustrative purposes.
xd(:,1) =3D sqrt(xd(:,1));
**************************************************************************
***************************** nhistmax.m
*********************************
function [nhist, m0] =3D nhistmax(x, nbins, low, high)
%NHISTMAX Compute normalized histogram.
% This function computes normalized unbiased histogram
% of x distribution and finds maximum peak magnitude.
% Histogram of x distribution among nbins
hx =3D hist(x,1:nbins);
% Compute normalized unbiased histogram
unbias =3D nbins:-1:1;
nhist =3D hx./unbias;
% Find maximal value between low and high bounds
m0 =3D max(nhist(low:high));
**************************************************************************
*************************** histsmooth.m
*********************************
function xs =3D histsmooth(x)
%HISTSMOOTH Smooth a histogram.
% Simple 3-point moving-average smoothing:
% xs(i) =3D a1*x(i-1) + a0*x(i) + a1*x(i+1)
%
% (This function can be replaced with any suitable
% histogram smoothing procedure)
a0 =3D 0.5;
a1 =3D 0.25;
nx=3Dlength(x);
xl =3D [x(2:nx) x(nx)];
xr =3D [x(1) x(1:(nx-1))];
xs =3Da0*x + a1*(xl+xr);
**************************************************************************
************************* localpeaks.m
***********************************
function peaks =3D localpeaks(hx, low, high, minsep)
%LOCALPEAKS Find positions and magnitudes of local peaks.
% This function searches for all local maxima in
% a histogram between low and high search bounds.
% INPUT:
% hx - input histogram vector
% low - lower bound to search for peaks
% high - upper bound to search for peaks
% minsep - minimal allowed separation between adjacent peaks
%
% OUTPUT:
% peaks =3D [pv pi] - Output array with 2 columns:
% pv has peak magnitudes from highest to lowest,
% pi has corresponding peak positions (in samples).
nx=3Dlength(hx); % total number of bins in a histogram
d=3Dsign(diff(hx));
peaks =3D [];
flag =3D 0;
% Find all local maxima between low and high bounds
for i=3Dlow:high-1
if (d(i) =3D=3D 1.0)
flag =3D 1;
nstart =3D i;
elseif (d(i) =3D=3D -1.0)
if (flag =3D=3D 1)
flag =3D 0;
nend =3D i;
pi =3D i - floor((nend-nstart)/2.0);
pv =3D hx(pi);
peaks =3D [peaks; pv pi];
end
end
end
[npeaks, junk] =3D size(peaks); % total number of peaks found
if (npeaks > 0)
peaks=3Dflipud(sortrows(peaks)); % Sort based on peak heights
% Optionally scan all peaks from highest to lowest and retain only
% those separated by at least minsep samples from already selected.
if (minsep > 2)
seppeaks =3D [peaks(1,:)]; % select highest peak first
for i=3D2:npeaks % analyze remaining peaks
[np, junk] =3D size(seppeaks);
count=3D0;
for j=3D1:np
if(abs(peaks(i,2)-seppeaks(j,2)) < minsep)
count =3D count+1;
end
end
if (count =3D=3D 0)
seppeaks =3D [seppeaks; peaks(i,:)];
end
end
peaks =3D seppeaks;
end
end
**************************************************************************
************************** svdembed.m
************************************
function xe =3D svdembed(x,ndim,nsvd)
%SVDEMBED Embed a signal using SVD-embedding procedure.
% This is an alternative to time-delay embedding.
%
% x - input signal vector
% ndim - embedding dimension
% nsvd - SVD window size in samples
% xe - output matrix with ndim columns and nr rows
% (Each column represents a separate dimension and
% each row represents a vector in state space)
nx =3D length(x);
av=3Dmean(x);
x=3Dx-av;
nr =3D nx - (nsvd-1);
xs=3Dzeros(nr,nsvd);
for i =3D 1:nsvd
xs(:,i) =3D x(i:(nr+i-1));
end
[U,S,V] =3D svd(xs,0);
Vd =3D V(:,1:ndim);
xe =3D xs*Vd;
[mr, junk] =3D size(xe);
% Normalize each column so that each dimension is in [0 1] range
for i =3D 1:ndim
xmax =3D max(xe(:,i));
xmin =3D min(xe(:,i));
xrange =3D xmax-xmin;
if (xrange > 0)
xe(:,i) =3D (xe(:,i) - xmin)/xrange;
else
xe(:,i) =3D zeros(mr,1);
end
end
**************************************************************************
robert bristow-johnson wrote:
> fizteh89 wrote:
> > Robert, not to offend you or anybody else in this group: by "they" I
> > certainly didn't mean some regulars here or some mom-and-pop workshops
> > out there.
> > "They" know who they are and what they are doing all too well...
>
> impressive.
>
> > I just downloaded and unzipped the code and it works fine for me, aside
> > from giving some warning messages. You need Windows version of Matlab
>
> i wrote:
> > > i downloaded pdemo and unzipped to a PC with MATLAB, changed the
> > > directory, and typed in pdemo at the prompt. ...
>
> > though and mine is pretty old, but it should be backward compatible.
>
> and MATLAB preaches strict backward compatibility. why does my windows
> version of MATLAB say the following?:
>
> > > =BB pdemo
> > > ??? c:\program files\matlab\pdemo\pdemo.p is not a MATLAB P-file.
>
> i'm not averse to fixing something in my MATLAB to get it to work, but
> i have no idea why it would do that. i download, unzip, and run it
> just as you web page says, and i can't get it to work.
>
> > I am sure other people out there were able to run the demo. "They" also
> > know how to search IFW of a pending US patent application to get all
> > the public info, including the source code submitted with provisional
> > application...
> > BTW, I posted some Matlab source here, in comp.dsp, some time ago.
>
> i started at
> http://groups.google.com/groups/profile?enc_user=3DriwiGBQAAAC5hnaoSNtzWu=

2Wz9LFDnOUOPANdqfI6prRsqjc7uCt1A
>
> and couldn't find it. i found some interesting quotes:
>
> Re: 23k$US for Matlab Embedded DSP development?
> fizteh89 wrote:
> > The world is not flat as far as local girls and local restaurants are
> > concerned.
> >
> > But, of course, some people prefer cyberfood or cybersex...
>
>
> > As for your "unanswered questions"...,well, you promised to come up
> > with some waveform to fool my PDA, but the truth is: there is no such
> > waveform.
>
> i am not going to go through the effort of implementing the algorithm
> that is described in your paper. if i can get a simple script from you
> that does only what is fully described in your paper (that leaves the
> SVD out, at least until the next pass) and not including any trade
> secreted stuff that might be in your pdemo.p file if it worked.
>
> if you do that, i will test it with waveforms i can either find or
> generate and give you (and post) an answer of how it went. if i find
> a waveform that messes it up, i will post the MATLAB script that made
> it, or if it is a recorded pitched signal, i'll email it to anyone who
> asks (in reasonable quantity).
>
> > ANY (benigh) periodic waveform you can come up with will produce
> > perfect pitch peaks - at 1, 2, 3... pitch periods - an "ideal" pitch
> > detector indeed.
>
> what if the peak is a little better at 2 or 3 than at one? which one
> do you pick? always the one with the maximum histogram? if that is
> the case, i know exactly how to construct a waveform that you think is
> at period 1 if you listened to it, but will show peak 2 (at twice the
> lag as peak 1) to be the maximum histogram (or the maxium of any other
> algorithm that picks the true and smallest period of a perfectly period
> signal in zero noise). to choose peak 1 over peak 2 when an inuadible
> but non-zero subharmonic (say one octave below the pitch you think you
> hear) requires something else that is not in your paper. if you always
> pick the very best peak in that histogram in your paper, it is
> susceptable to the same kind of octave error problem other PDAs have
> (unless there's something else there other than picking the "true" or
> "mathematically best" period).
>
> > This is as good as it gets - there will never be anything better for
> > pitch/fundamental period detection.
>
> you should sell cars. you're really good at it.
>
> > And your favorite example of mixing very weak 100 Hz sine wave with
> > strong 200 Hz sine wave just proves the validity of the method: the
> > true fundamental frequency is 100 Hz and the method is more sensitive
> > than human ear, subject only to a graceful degradation of accuracy with
> > increasing signal-to-noise ratio...
>
> whoa! woo-haa! that'll teach me not to read to the bottom of the post
> before replying! (actually it hadn't in the past.)
>
> so you're saying that if you were listening to a pure sine wave
> synthesized by the pitch information that comes from your pitch
> detector as it is tracking some given input like a human voice in
> speech or singing, and if an inaudible subharmonic (say -80 dB of the
> main tone) somehow finds its way into your tone, the PDA (and your
> synthesized sine wave) will immediately jump down an octave and then
> toggle back up when the subharmonic is gone (possibly multiple times
> during the tone duration)? if that is the case, i would say that your
> PDA is not useful for tough pitch detection applications. that
> behavior represents a failure of correct operation of a PDA. it is an
> "octave error", a glitch of some manner.
>
> > I hope that I don't sound paternalistic or arrogant.
>
> no, of course not.
>
> > After all, it's
> > just a very small (yet essential) part of DSP and I DID spent a LOT of
> > time and effort to work it out in sufficient detail for DSP practitione=

rs.
>
> that you have spent a lot of time working it out (and writing the paper
> and filing the patent and creating the website and posting the pdemo.p
> file that i can't get to work in MATLAB) is not the issue. nobody is
> disputing that.
>=20
> r b-j