kpitch.cpp

这个是用来进行基频检测的

// KPitch.cpp: implementation of the KPitch class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "KPitch.h"
#include <math.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

KPitch::KPitch():m_win(NULL),m_lift(NULL)
{
}

KPitch::~KPitch()
{
	Destroy();
}

void KPitch::Destroy()
{
	if(m_win){ delete [] m_win; m_win = NULL; }
	if(m_lift){ delete [] m_lift; m_lift = NULL; }
	if(m_buff1) { delete [] m_buff1; m_buff1 = NULL; } 
	if(m_buff2) { delete [] m_buff2; m_buff2 = NULL; }
	if(m_buff3) { delete [] m_buff3; m_buff3 = NULL; }
	if(m_buff4) { delete [] m_buff4; m_buff4 = NULL; }
}

void KPitch::Create(int fs, int fmin, int fmax)
{
	int i;
	m_fs = fs;
	m_nmin = fs/fmax;
	m_nmax = fs/fmin;
	m_winsize = int(fs*3/fmin);
	m_fftsize = int(pow(2,int(log(m_winsize*2-1)/log(2)))*2);
	m_win = new double[m_winsize];
	for(i=0;i<m_winsize;i++) // hanning window
		m_win[i] = 0.5 - 0.5 * cos(2*pi*(i+0.5)/m_winsize);
	m_lift = new double[m_fftsize];
	m_buff1 = new double[m_fftsize];
	m_buff2 = new double[m_fftsize];
	m_buff3 = new double[m_fftsize];
	m_buff4 = new double[m_fftsize];
	for(i=0;i<m_winsize;i++)
		m_lift[i] = m_win[i];
	for(i=m_winsize;i<m_fftsize;i++)
		m_lift[i] = 0;
	for(i=0;i<m_fftsize;i++)
		m_buff1[i] = 0;
	fft(m_lift,m_buff1,m_buff2,m_buff3,m_fftsize);
	for(i=0;i<m_fftsize;i++)
		m_buff2[i] = m_buff2[i]*m_buff2[i]+m_buff3[i]*m_buff3[i];
	for(i=0;i<m_fftsize;i++)
		m_buff3[i] = 0;
	fft(m_buff2,m_buff3,m_lift,m_buff1,m_fftsize);
	for(i=1;i<m_fftsize;i++)
		m_lift[i] /= m_lift[0];
	m_lift[0] = 1.0;
	for(i=0;i<m_fftsize;i++)
		m_lift[i] = sqrt(m_lift[i]);
	for(i=0;i<m_fftsize;i++)
		if(m_lift[i]!=0.0)
			m_lift[i] = 1.0/m_lift[i];
}

void KPitch::fft(double rin[], double iin[], double rout[], double iout[], int n)
{
	/*
	   n must be power of 2
	*/
	int i,j,s,t,m,offset,logn;
	double wr,wi,qr,qi;
	
	// bit invert
	logn = int(log(n+0.5)/log(2));
	for(i=0;i<n;i++)
	{
		s=i; t=0; for(j=0;j<logn;j++){ t*=2; t+=s&1; s/=2;}
		rout[i] = rin[t];
		iout[i] = iin[t];
	}
	// base 2-->4-->8-->..-->n
	m=2;
	while(m<=n)
	{
		s=n/m;
		for(i=0;i<s;i++)
		{
			offset = i*m;
			for(j=0;j<m/2;j++)
			{
				wr = cos(2*pi*j/m);
				wi = sin(2*pi*j/m);

				qr = rout[offset+j+m/2]*wr-iout[offset+j+m/2]*wi;
				qi = iout[offset+j+m/2]*wr+rout[offset+j+m/2]*wi;

				rout[offset+j+m/2] = rout[offset+j] - qr;
				iout[offset+j+m/2] = iout[offset+j] - qi;
				rout[offset+j] += qr;
				iout[offset+j] += qi;
			}
		}
		m*=2;
	}
}

void KPitch::Wav2Pitch(double wav[], int nwav, double hoptime, double f0[], int nf0)
{
	double maxv;
	int i,j,p,offset,frmcnt;
	frmcnt = GetFrameCount(nwav,hoptime);
	if(frmcnt<nf0)
	{
		for(i=frmcnt;i<nf0;i++)
			f0[i] = 0.0;
		nf0 = frmcnt;
	}
	for(i=0;i<nf0;i++)
	{
		offset = int(hoptime*i);
		// get correlation in m_buff1[j]
		for(j=0;j<m_winsize;j++)
			m_buff1[j] = wav[j+offset]*m_win[j];
		for(;j<m_fftsize;j++)
			m_buff1[j] = 0.0;
		for(j=0;j<m_fftsize;j++)
			m_buff2[j] = 0.0;
		fft(m_buff1,m_buff2,m_buff3,m_buff4,m_fftsize);
		for(j=0;j<m_fftsize;j++)
			m_buff3[j] = m_buff3[j]*m_buff3[j]+m_buff4[j]*m_buff4[j];
		for(j=0;j<m_fftsize;j++)
			m_buff4[j] = 0.0;
		fft(m_buff3,m_buff4,m_buff1,m_buff2,m_fftsize);
		if(m_buff1[j]!=0.0)
		{
			for(j=1;j<m_fftsize;j++)
				m_buff1[j] /= m_buff1[0];
			m_buff1[0] = 1.0;
			for(j=0;j<m_fftsize;j++)
				m_buff1[j] *= m_lift[j];
			// find max peak
			p=0;maxv=0;
			for(j=m_nmin;j<m_nmax;j++)
				if( (m_buff1[j]>maxv) && 
					(m_buff1[j]>=m_buff1[j-1]) && 
					(m_buff1[j]>=m_buff1[j+1]))
				{ p = j; maxv = m_buff1[j]; }
			// get pitch
			if(p>0)
				f0[i] = m_fs*1.0/p;
			else
				f0[i] = 0.0;
		}else{
			f0[i] = 0.0;
		}
	}
}

int KPitch::GetFrameCount(int nwav, double hoptime)
{
	return (nwav-int(hoptime*m_fs))/m_winsize + 1;
}