Hi, I want to show a Pitch graph/curve of a wave file. I saw that good results can be made by using AutoCorrelation. I have the autocorrelation function and I also have the FFT function, but I don't know what I should put into the autocorrelation function in order to get a pitch line. I have a wave file (spoken text), and also have the FFT results of the wave file, how can I extract from this data a nice Pitch curve. The purpose of the pitch line is to show the intonations of the text. If anybody could explain to me the steps to be taken in plain english, pseudo code, or C#, VB(.NET) or Java, I would be very very thankfull. Many Thanks, Gidon
Pitch calculation (FFT or autocorrelation)
Started by ●June 5, 2006
Reply by ●June 5, 20062006-06-05
gjunge wrote:> Hi, > > I want to show a Pitch graph/curve of a wave file. I saw that good results > can be made by using AutoCorrelation. I have the autocorrelation function > and I also have the FFT function, but I don't know what I should put into > the autocorrelation function in order to get a pitch line. > I have a wave file (spoken text), and also have the FFT results of the > wave file, how can I extract from this data a nice Pitch curve. The > purpose of the pitch line is to show the intonations of the text. > > If anybody could explain to me the steps to be taken in plain english, > pseudo code, or C#, VB(.NET) or Java, I would be very very thankfull. > > Many Thanks, > GidonI am not sure what you are asking. Do you want pitch as a function of time?
Reply by ●June 6, 20062006-06-06
Neither FFT nor autocorrelation will do nowadays... Pitch detection problem was laid to rest a few years ago... Check this page: http://www.soundmathtech.com/pitch And read US Patent Application 20030088401 at www.uspto.gov/patft (Strictly non-commercial license applies) That's what they use for pitch (fundamental frequency) detection nowadays, but don't ask - nobody will admit it :-) gjunge wrote:> Hi, > > I want to show a Pitch graph/curve of a wave file. I saw that good results > can be made by using AutoCorrelation. I have the autocorrelation function > and I also have the FFT function, but I don't know what I should put into > the autocorrelation function in order to get a pitch line. > I have a wave file (spoken text), and also have the FFT results of the > wave file, how can I extract from this data a nice Pitch curve. The > purpose of the pitch line is to show the intonations of the text. > > If anybody could explain to me the steps to be taken in plain english, > pseudo code, or C#, VB(.NET) or Java, I would be very very thankfull. > > Many Thanks, > Gidon
Reply by ●June 6, 20062006-06-06
fizteh89 wrote:> Neither FFT nor autocorrelation will do nowadays... > > Pitch detection problem was laid to rest a few years ago... > > Check this page: http://www.soundmathtech.com/pitchthis website says: "Currently the demonstration software is provided in the form of several Matlab P-files. Download the compressed ZIP-archive pdemo.zip (151 Kb) Unpack it in the Matlab work directory. Then type pdemo at the Matlab command prompt. You should see the GUI window coming up. You can now load and process speech files in wav or ascii formats (Clean speech sampled at 16 kHz is assumed)." i downloaded pdemo and unzipped to a PC with MATLAB, changed the directory, and typed in pdemo at the prompt. this is what i get: =BB pwd ans =3D C:\Program Files\MATLAB\pdemo =BB dir . nhist.p selpeaks.p xdistall.p .=2E pdemo.p shehadyourdark.wav embed.p pframe.p showframe.p histsm.p ptrack.p xdist.p =BB pdemo ??? c:\program files\matlab\pdemo\pdemo.p is not a MATLAB P-file. is there something i'm doing wrong? it doesn't work. i've read your paper, you have read my response to it years ago: http://groups.google.com/group/comp.dsp/browse_thread/thread/5f5b466f059eb9= ce http://groups.google.com/group/comp.dsp/msg/9daa9e1f9a0dcba0 and i have yet to hear how well this pitch detector works and sounds.> And read US Patent Application 20030088401 at www.uspto.gov/patft > (Strictly non-commercial license applies) > > That's what they use for pitch (fundamental frequency) detection > nowadays, but don't ask - nobody will admit it :-)who is "they"? if there are all sorts of unidentified people or companies using your algorithm that will not say so, what meaning can we ascribe to the statement above? Dmitry, i am not writing your PDA off, i just haven't been able to see anything concrete other than your paper http://www.soundmathtech.com/pitch/paper.html and the issues i have brought up regarding that have not been addressed. best regards, r b-j
Reply by ●June 6, 20062006-06-06
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... 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 though and mine is pretty old, but it should be backward compatible. 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. 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. 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. This is as good as it gets - there will never be anything better for pitch/fundamental period detection. 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... I hope that I don't sound paternalistic or arrogant. 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 practitioners. robert bristow-johnson wrote:> fizteh89 wrote: > > Neither FFT nor autocorrelation will do nowadays... > > > > Pitch detection problem was laid to rest a few years ago... > > > > Check this page: http://www.soundmathtech.com/pitch > > this website says: > > "Currently the demonstration software is provided in the form of > several Matlab P-files. > > Download the compressed ZIP-archive pdemo.zip (151 Kb) > > Unpack it in the Matlab work directory. Then type pdemo at the Matlab > command prompt. You should see the GUI window coming up. You can now > load and process speech files in wav or ascii formats (Clean speech > sampled at 16 kHz is assumed)." > > > > > i downloaded pdemo and unzipped to a PC with MATLAB, changed the > directory, and typed in pdemo at the prompt. this is what i get: > > =BB pwd > > ans =3D > > C:\Program Files\MATLAB\pdemo > > =BB dir > > . nhist.p selpeaks.p xdistall.p > > .. pdemo.p shehadyourdark.wav > embed.p pframe.p showframe.p > histsm.p ptrack.p xdist.p > > =BB pdemo > ??? c:\program files\matlab\pdemo\pdemo.p is not a MATLAB P-file. > > is there something i'm doing wrong? it doesn't work. > > i've read your paper, you have read my response to it years ago: > > http://groups.google.com/group/comp.dsp/browse_thread/thread/5f5b466f059e=b9ce> http://groups.google.com/group/comp.dsp/msg/9daa9e1f9a0dcba0 > > and i have yet to hear how well this pitch detector works and sounds. > > > And read US Patent Application 20030088401 at www.uspto.gov/patft > > (Strictly non-commercial license applies) > > > > That's what they use for pitch (fundamental frequency) detection > > nowadays, but don't ask - nobody will admit it :-) > > who is "they"? if there are all sorts of unidentified people or > companies using your algorithm that will not say so, what meaning can > we ascribe to the statement above? > > Dmitry, i am not writing your PDA off, i just haven't been able to see > anything concrete other than your paper > http://www.soundmathtech.com/pitch/paper.html and the issues i have > brought up regarding that have not been addressed. >=20 > best regards, >=20 > r b-j
Reply by ●June 6, 20062006-06-06
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 Matlabi 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=3DriwiGBQAAAC5hnaoSNtzWu2W= z9LFDnOUOPANdqfI6prRsqjc7uCt1A 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 practitioners.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. r b-j
Reply by ●June 7, 20062006-06-07
gjunge wrote:> I want to show a Pitch graph/curve of a wave file. I saw that good results > can be made by using AutoCorrelation. I have the autocorrelation function > and I also have the FFT function, but I don't know what I should put into > the autocorrelation function in order to get a pitch line. > I have a wave file (spoken text), and also have the FFT results of the > wave file, how can I extract from this data a nice Pitch curve. The > purpose of the pitch line is to show the intonations of the text.Google for "phase vocoder" for how to use FFT results of a series of windows of data offset in time; and for "cepstral" or "cepstrum" methods to help with pitch determination of ambiguous sets of frequencies. IMHO. YMMV. -- Ron N. http://www.nicholson.com/rhn
Reply by ●June 7, 20062006-06-07
#include "stdafx.h"
#define W_max 256 // Max. Number of window length
#define PI 3.14159265358979
#define W_length 256 // window length
#define S_length 128 // shift length
#define P 256 // r(k) order
// Function headers
void CalcAuto(double *pX, double *pRxx);
void FrameProcessing(int index, short *pData, double *pHcoef, double
*pX);
double GetPitch(double *pRxx, int Fs);
// Global variable
FILE *stream;
void main(void)
{
CWaveSoundRead *pWaveFileReader;
// =EC=9B=A8=EC=9D=B4=EB=B8=8C =ED=8C=8C=EC=9D=BC =EA=B4=80=EB=A0=A8 =EB=
=B3=80=EC=88=98
int nFileSize =3D 0; // wav file=EC=9D=98 =ED=81=AC=EA=B8=B0
int SamplingRate =3D 0; // wav file=EC=9D=98 sampling rate
int Resolution =3D 0; // wav file=EC=9D=98 resolution
short *BUFFER =3D NULL; // =EB=B2=84=ED=8D=BC (=EC=9B=A8=EC=9D=B4=
=EB=B8=8C =EB=8D=B0=EC=9D=B4=ED=84=B0
=EC=A0=80=EC=9E=A5=EC=9A=A9)
pWaveFileReader =3D NULL;
BYTE *pbWavData; // Pointer to actual wav data
UINT cbWavSize; // Size of data
char fileName[20]; // =EC=9B=A8=EC=9D=B4=EB=B8=8C =ED=8C=8C=EC=
=9D=BC =EC=9D=BD=EA=B8=B0=EC=9A=A9
=ED=8C=8C=EC=9D=BC=EC=9D=B4=EB=A6=84
// Pitch extraction =EA=B4=80=EB=A0=A8 =EB=B3=80=EC=88=98
int tframe =3D 0; // Number of frames
int datalength =3D 0; // data length.
double pitch =3D 0; // pitch =EA=B3=84=EC=82=B0=EA=B2=B0=EA=B3=
=BC
double h[W_length]; // Hamming window coef.
double x[W_length]; // Hamming windowed data
double Rxx[W_length]; // autocorrelation
double *Pitch; // fundamental freq. (pitch)
sprintf(fileName, "communication.wav");
// =EA=B3=84=EC=82=B0=EA=B0=92 =ED=99=95=EC=9D=B8=EC=9C=84=ED=95=9C file p=
rocessing
stream =3D fopen("Result.txt", "w");
pWaveFileReader=3Dnew CWaveSoundRead();
// Wav file =EC=97=B4=EA=B8=B0.
if(FAILED( pWaveFileReader->Open((LPSTR)(LPCTSTR)fileName) ))
{
cout<<"Failed to Open the wave file"<<endl;
exit(0);
}
nFileSize =3D pWaveFileReader->m_ckIn.cksize; // =ED=8C=8C=EC=9D=BC =
=ED=81=AC=EA=B8=B0
pbWavData =3D new BYTE[nFileSize]; // BYTE: 8 bits
BUFFER =3D new short[nFileSize/2]; // =EB=B2=84=ED=8D=BC
// =EC=9B=A8=EC=9D=B4=EB=B8=8C =ED=8C=8C=EC=9D=BC =EC=9D=BD=EA=B8=B0
if(FAILED( pWaveFileReader->Read(nFileSize, pbWavData, &cbWavSize) ))
{
cout<<"Failed to Read the wave file"<<endl;
exit(0);
}
pWaveFileReader->Reset(); // file=EC=9D=84 =EC=B4=88=EA=B8=B0 =EC=
=9C=84=EC=B9=98=EB=A1=9C
=EC=B4=88=EA=B8=B0=ED=99=94.
SamplingRate =3D pWaveFileReader->m_pwfx->nSamplesPerSec; // wav file=EC=
=9D=98
sampling rate.
Resolution =3D pWaveFileReader->m_pwfx->wBitsPerSample; // wav file=EC=9D=
=98
resolution.
// =EB=B2=84=ED=8D=BC=EC=97=90 =EC=A0=80=EC=9E=A5
memcpy(BUFFER, pbWavData, nFileSize);
// wav file =EC=A0=95=EB=B3=B4 =EC=B6=9C=EB=A0=A5
fprintf(stream, "Sampling rate: %d \n", SamplingRate );
fprintf(stream, "Resoultion : %d \n", Resolution );
fprintf(stream, "File size : %d \n", nFileSize );
fprintf(stream, "Size of data : %d \n", cbWavSize );
fprintf(stream,
"--------------------------------------------------------- \n" );
// frame =EA=B0=AF=EC=88=98 =EA=B3=84=EC=82=B0
datalength =3D nFileSize/2; // wav file=EC=9D=98 =EB=8D=B0=EC=9D=
=B4=ED=84=B0 =EA=B0=AF=EC=88=98
tframe =3D (int)((datalength-W_length)/S_length) + 1;
// =EA=B0=81 frame=EC=9D=98 pitch=EB=A5=BC =EC=A0=80=EC=9E=A5=ED=95=98=EA=
=B8=B0=EC=9C=84=ED=95=9C =EB=B3=80=EC=88=98
Pitch =3D new double[tframe];
for (int i=3D0; i<W_length; i++)
{
h[i] =3D 0.54 - 0.46*cos(2*PI*i/(double)(W_length-1));
fprintf(stream, "h[%d]: %f \n", i, h[i]);
}
// =EB=A7=A4 frame =EB=A7=88=EB=8B=A4 pitch =EA=B2=80=EC=B6=9C=ED=95=A8.
for (int frame=3D0; frame<tframe; frame++)
{
// Frame processing
FrameProcessing(frame, BUFFER, h, x);
// Autocorrelation analysis
CalcAuto(x, Rxx);
// Pitch extraction
Pitch[frame] =3D GetPitch(Rxx, SamplingRate);
}
delete pbWavData;
delete pWaveFileReader;
pWaveFileReader->Close();
pWaveFileReader=3DNULL;
delete [] Pitch;
fclose(stream);
}
/*--------------------------------------------------------------------
double GetPitch(double *pRxx, // autocorrelation
int Fs) // sampling rate
--------------------------------------------------------------------*/
double GetPitch(double *pRxx,int Fs)
{
double s1 =3D 0; // slope 1
double s2 =3D 0; // slope 2
double peakValue =3D 0; // peak=EC=A0=90=EC=97=90=EC=84=9C=EC=9D=
=98 r(k) =EA=B0=92.
double peakIndex =3D 0; // peak index.
double energy =3D 0; // short-time energy (r(0))
double f0_temp =3D 0; // temporary f0
double f0 =3D 0; // fundamental freq. (pitch)
energy =3D 10 * log10(pRxx[0]);
// speech =EA=B5=AC=EA=B0=84=EC=9D=BC =EA=B2=BD=EC=9A=B0.
if (energy > 40)
{
// pre-emphasis factor =EA=B3=84=EC=82=B0. voiced/unvoiced =EA=B5=AC=EB=
=B6=84
double mu =3D pRxx[1]/pRxx[0];
// voiced=EC=9D=BC =EA=B2=BD=EC=9A=B0.
if (mu > 0.5)
{
// peak-picking
for (int k=3D0; k<W_length-1; k++)
{
s1 =3D pRxx[k+1] - pRxx[k];
s2 =3D pRxx[k+2] - pRxx[k+1];
// =EA=B8=B0=EC=9A=B8=EA=B8=B0=EB=A5=BC =EC=9D=B4=EC=9A=A9 =EC=B5=9C=EB=
=8C=80 peak =EB=B0=8F peak index =EA=B2=80=EC=B6=9C.
if (s1>0 && s2<0)
{
if (pRxx[k+1] > peakValue)
{
peakValue =3D pRxx[k+1];
peakIndex =3D k+1;
}
}
}
// pitch =EA=B3=84=EC=82=B0
if (peakIndex-1 > 0)
{
f0_temp =3D Fs/(peakIndex-1);
}
// =EB=82=98=EB=88=97=EC=85=88=EC=9D=B4 0=EC=9D=B4 =EB=90=98=EB=8A=94 =
=EC=A1=B0=EA=B1=B4.
else if (peakIndex-1 =3D=3D 0)
{
f0_temp =3D 0;
}
// pitch =EC=A1=B4=EC=9E=AC=EB=B2=94=EC=9C=84=EC=97=90 =EC=9E=88=EC=9D=
=84 =EA=B2=BD=EC=9A=B0=EC=97=90=EB=A7=8C f0 =EC=84=A4=EC=A0=95 (40~400 Hz)
if (f0_temp < 40)
{
f0 =3D 0;
}
else if (f0_temp > 400)
{
f0 =3D 0;
}
else
{
f0 =3D f0_temp;
}
}
// unvoiced=EC=9D=BC =EA=B2=BD=EC=9A=B0.
else
{
f0 =3D 0;
}
}
// silence =EA=B5=AC=EA=B0=84=EC=9D=BC =EA=B2=BD=EC=9A=B0.
else
{
f0 =3D 0;
}
return f0;
}
/*--------------------------------------------------------------------
void CalcAuto(double *pX, // windowed data
double *pRxx) // autocorrelation
--------------------------------------------------------------------*/
void CalcAuto(double *pX, double *pRxx)
{
int k; // r(k)=EC=9D=98 index.
int n;
for (k=3D0; k<P; k++)
{
pRxx[k] =3D 0.0;
for (n=3D0; n<W_length-k; n++)
{
pRxx[k] =3D pRxx[k] + pX[n]*pX[n+k];
}
pRxx[k] =3D pRxx[k]/W_length;
}
}
/*--------------------------------------------------------------------
void FrameProcessing(int index, // Frame index
short *pData, // Wav data
double *pHcoef, // Hamming window coef.
double *pX) // Hamming windowed data
--------------------------------------------------------------------*/
void FrameProcessing(int index, short *pData, double *pHcoef, double
*pX)
{
// Hamming Windowing & rectangular windowing
/*
pX[0] =3D pHcoef[0] * (double)pData[index*S_length];
for (int i=3D1; i<W_length; i++)
{
pX[i] =3D pHcoef[i] * (double)pData[index*S_length + i];
}
*/
// Rectangular windowing
pX[0] =3D (double)pData[index*S_length];
for (int i=3D1; i<W_length; i++)
{
pX[i] =3D (double)pData[index*S_length + i];
}
}
gjunge wrote:
> Hi,
>
> I want to show a Pitch graph/curve of a wave file. I saw that good results
> can be made by using AutoCorrelation. I have the autocorrelation function
> and I also have the FFT function, but I don't know what I should put into
> the autocorrelation function in order to get a pitch line.
> I have a wave file (spoken text), and also have the FFT results of the
> wave file, how can I extract from this data a nice Pitch curve. The
> purpose of the pitch line is to show the intonations of the text.
>
> If anybody could explain to me the steps to be taken in plain english,
> pseudo code, or C#, VB(.NET) or Java, I would be very very thankfull.
>=20
> Many Thanks,=20
> Gidon
Reply by ●June 7, 20062006-06-07
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
Reply by ●June 8, 20062006-06-08
fizteh89 wrote:> OK, Robert, just for you... I am posting it a second time.i missed it the first time. and i couldn't find it with Google groups, but i didn't go through each of your posts.> This is cut-and-paste Matlab routine - just fix the wrapped around > comment lines.i can do that. (and i will within a week i think.)> Add some signal pre-processing (like low-pass filtering) as you wish...i can do that, too. most likely i will couple your PDA output to an oscillator and synthesize the note that is tracking the input. maybe your code does it, i haven't examined it yet.> You are still confused, my friend, about pitch vs. fundamental > frequency (or, rather, fundamental period)i think i know the meaning of all three.> Pitch is an auditory concept, it's tied to human perception.no argument from me about that.> 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.i would define "PDA" to be closer to what the acronym stands for, a "Pitch Detection Algorithm". it is not soley a f0 estimator although there is a very strong relationship between f0 and pitch. the problem is that with many instruments and the human voice, both spoken and singing, sometimes the nature of the vibration of whatever it is (reed, lips, strings, vocal cords, or a column of air) is that one vibration is nearly identical to the next but slightly different but the following vibration is very much like the first. every odd cycle is more identical to each other than they are to the even cycles. so the f0 is truly one octave lower, but the human listening may not perceive it to be that lower octave if the even cycles are soooo similar to the odd cycles that the difference is hard to audibly measure (but a machine could, like a PDA). anyway, any application i can think of for a PDA desires the pitch returned to be the same as the perceived pitch. if the note sounds like middle C, the PDA should say it's middle C, not the C an octave below nor the F that is 7 semitones lower than that. if the PDA has trouble making up its mind about whether it is the low C or middle C, that will sound like crap when the PDA output is controlling an NCO (numerically controlled oscillator) because the NCO will jump back and forth between the different notes spuriously.> (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 informationwhat does the following sound like?:>>> This is as good as it gets - there will never be anything better for >>> pitch/fundamental period detection.... a scientist or a salesman? r b-j
