#include "waveconvertor.h"
//
// class CWaveConvertor
/* Pjotr '87. */
// As best as I can determine this is in the public domain. No license was found
// withit what so ever. The author is listed as:
// Peter Valkenburg (valke@cs.vu.nl).
//
// If there is some disagreement then contact Bruce Forsberg (forsberg@tns.net)
#include
#define pi 3.1415926535897932384626434
#define c_re(c) ((c).re)
#define c_im(c) ((c).im)
/* C_add_mul adds product of c1 and c2 to c. */
#define c_add_mul(c, c1, c2) { COMPLEX C1, C2; C1 = (c1); C2 = (c2); \
c_re (c) += C1.re * C2.re - C1.im * C2.im; \
c_im (c) += C1.re * C2.im + C1.im * C2.re; }
/* C_conj substitutes c by its complex conjugate. */
#define c_conj(c) { c_im (c) = -c_im (c); }
/* C_realdiv divides complex c by real. */
#define c_realdiv(c, real) { c_re (c) /= (real); c_im (c) /= (real); }
/*
* W gives the (already computed) Wn ^ k (= e ^ (2pi * i * k / n)).
* Notice that the powerseries of Wn has period Nfactors.
*/
#define W(n, k) (W_factors [((k) * (Nfactors / (n))) % Nfactors])
/*! \brief Constructor.
*/
aflibFFT::aflibFFT()
{
Nfactors = 0;
W_factors = NULL;
}
/*! \brief Destructor.
*/
aflibFFT::~aflibFFT()
{
if (W_factors != NULL)
delete W_factors;
}
/*! \brief Performs a forward or reverse FFT.
This is the main API is this class. It will perform either a forward or
inverse FFT depending how InverseTransform is set. If set to FALSE then
forward FFT will be performed, TRUE and a inverse FFT will be performed.
The number of samlpes (NumSamples) must be a power of 2. The ImagIn
pointer can be NULL if there are no imaginary values. The user is
responsable for passing in pointers for RealOut and ImagOut containing
arrays of the proper size.
*/
void
aflibFFT::fft_double (
unsigned NumSamples,
int InverseTransform,
const double *RealIn,
const double *ImagIn,
double *RealOut,
double *ImagOut )
{
COMPLEX in[1024];
COMPLEX out[1024];
COMPLEX * in_local = NULL;
COMPLEX * out_local = NULL;
register COMPLEX * in_ptr;
register COMPLEX * out_ptr;
register unsigned int i;
// IF 1024 samples or less use local buffer else allocate memory
if (NumSamples > 1024)
{
in_local = new COMPLEX[NumSamples];
out_local = new COMPLEX[NumSamples];
in_ptr = in_local;
out_ptr = out_local;
}
else
{
in_ptr = in;
out_ptr = out;
}
// Fill real and imaginary array
for (i = 0; i < NumSamples; i++)
{
c_re(in_ptr[i]) = RealIn[i];
if (ImagIn == NULL)
c_im(in_ptr[i]) = 0.0;
else
c_im(in_ptr[i]) = ImagIn[i];
}
// Perform transform
if (InverseTransform == TRUE)
{
rft(in_ptr, NumSamples, out_ptr);
}
else
{
fft(in_ptr, NumSamples, out_ptr);
}
// Fill real and imaginary array
for (i = 0; i < NumSamples; i++)
{
RealOut[i] = c_re(out_ptr[i]);
ImagOut[i] = c_im(out_ptr[i]);
}
// Free memory if local arrays were not used
if (in_local != NULL)
delete [] in_local;
if (out_local != NULL)
delete [] out_local;
}
/*
* Forward Fast Fourier Transform on the n samples of complex array in.
* The result is placed in out. The number of samples, n, is arbitrary.
* The W-factors are calculated in advance.
*/
int
aflibFFT::fft (
COMPLEX *in,
unsigned n,
COMPLEX *out)
{
unsigned i;
for (i = 0; i < n; i++)
c_conj (in [i]);
if (W_init (n) == -1)
return -1;
Fourier (in, n, out);
for (i = 0; i < n; i++) {
c_conj (out [i]);
c_realdiv (out [i], n);
}
return 0;
}
/*
* Reverse Fast Fourier Transform on the n complex samples of array in.
* The result is placed in out. The number of samples, n, is arbitrary.
* The W-factors are calculated in advance.
*/
int
aflibFFT::rft (
COMPLEX *in,
unsigned n,
COMPLEX *out)
{
if (W_init (n) == -1)
return -1;
Fourier (in, n, out);
return 0;
}
/*
* Recursive (reverse) complex fast Fourier transform on the n
* complex samples of array in, with the Cooley-Tukey method.
* The result is placed in out. The number of samples, n, is arbitrary.
* The algorithm costs O (n * (r1 + .. + rk)), where k is the number
* of factors in the prime-decomposition of n (also the maximum
* depth of the recursion), and ri is the i-th primefactor.
*/
void
aflibFFT::Fourier (
COMPLEX *in,
unsigned n,
COMPLEX *out)
{
unsigned r;
if ((r = radix (n)) < n)
split (in, r, n / r, out);
join (in, n / r, n, out);
}
/*
* Give smallest possible radix for n samples.
* Determines (in a rude way) the smallest primefactor of n.
*/
unsigned
aflibFFT::radix (unsigned n)
{
unsigned r;
if (n < 2)
return 1;
for (r = 2; r < n; r++)
if (n % r == 0)
break;
return r;
}
/*
* Split array in of r * m samples in r parts of each m samples,
* such that in [i] goes to out [(i % r) * m + (i / r)].
* Then call for each part of out Fourier, so the r recursively
* transformed parts will go back to in.
*/
void
aflibFFT::split (
register COMPLEX *in,
register unsigned r,
register unsigned m,
register COMPLEX *out)
{
register unsigned k, s, i, j;
for (k = 0, j = 0; k < r; k++)
for (s = 0, i = k; s < m; s++, i += r, j++)
out [j] = in [i];
for (k = 0; k < r; k++, out += m, in += m)
Fourier (out, m, in);
}
/*
* Sum the n / m parts of each m samples of in to n samples in out.
* r - 1
* Out [j] becomes sum in [j % m] * W (j * k). Here in is the k-th
* k = 0 k n k
* part of in (indices k * m ... (k + 1) * m - 1), and r is the radix.
* For k = 0, a complex multiplication with W (0) is avoided.
*/
void
aflibFFT::join (
register COMPLEX *in,
register unsigned m,
register unsigned n,
register COMPLEX *out)
{
register unsigned i, j, jk, s;
for (s = 0; s < m; s++)
for (j = s; j < n; j += m) {
out [j] = in [s];
for (i = s + m, jk = j; i < n; i += m, jk += j)
c_add_mul (out [j], in [i], W (n, jk));
}
}
int
aflibFFT::W_init(unsigned n)
{
unsigned k;
if (n == Nfactors)
return 0;
if (Nfactors != 0 && W_factors != 0)
delete [] W_factors;
if ((Nfactors = n) == 0)
return 0;
if ((W_factors = new COMPLEX[n]) == NULL)
return -1;
for (k = 0; k < n; k++) {
c_re (W_factors [k]) = cos (2 * pi * k / n);
c_im (W_factors [k]) = sin (2 * pi * k / n);
}
return 0;
}
//
// 构造函数
CWaveConvertor::CWaveConvertor(void)
{
}
//
// 析构函数
CWaveConvertor::~CWaveConvertor(void)
{
}
//
// 对输入数据进行短时快速FFT变换，输入数据为已经加窗的数据
void CWaveConvertor::ConverToFFT(
unsigned int nSamples, // 样本数量
unsigned int nShorts, // 短时点数
const double* pRealIn, // 输入的实部，不可为空
double* pRealOut, // 输出的实部
double* pImageOut // 输出的虚部
)
{
aflibFFT libFFT;
// 计算所需变换次数
unsigned int nCount = nSamples / nShorts;
// 数组偏移
unsigned int nOffSet = 0;
// 初始化输出数据
memset(pRealOut, 0, sizeof(double) * nSamples);
memset(pImageOut, 0, sizeof(double) * nSamples);
for (unsigned int i = 0; i < nCount; i++)
{
nOffSet = i * nShorts;
// 进行快速FFT
libFFT.fft_double(
nShorts,
FALSE,
&pRealIn[nOffSet],
NULL,
&pRealOut[nOffSet],
&pImageOut[nOffSet]
);
}
}
//
// 对输入数据进行Reverse Fourier转化，输入数据为已经加窗的数据
void CWaveConvertor::ConvertToRFT(
unsigned int nSamples, // 样本数量
unsigned int nShorts, // 短时点数
const double* pRealIn, // 输入的实部，不可为空
double* pRealOut, // 输出的实部
double* pImageOut // 输出的虚部
)
{
aflibFFT libFFT;
// 计算所需变换次数
unsigned int nCount = nSamples / nShorts;
// 数组偏移
unsigned int nOffSet = 0;
// 初始化输出数据
memset(pRealOut, 0, sizeof(double) * nSamples);
memset(pImageOut, 0, sizeof(double) * nSamples);
for (unsigned int i = 0; i < nCount; i++)
{
nOffSet = i * nShorts;
// 进行快速FFT
libFFT.fft_double(
nShorts,
TRUE,
&pRealIn[nOffSet],
NULL,
&pRealOut[nOffSet],
&pImageOut[nOffSet]
);
}
}
//
// 获取输入信号的功率谱
void CWaveConvertor::ConvertToPowerSpectral(
unsigned int nSamples, // 样本数量
unsigned int nShorts, // 短时点数
const double* pRealIn, // 输入信号
double* pDataOut // 输出的功率谱
)
{
double* pRealOut = NULL;
double* pImageOut = NULL;
// 输入数据是否有效
if (nSamples > 0 && pRealIn != NULL)
{
// 为FFT输出数据分配空间
pRealOut = new double[nSamples];
pImageOut = new double[nSamples];
// 进行FFT变换
CWaveConvertor::ConverToFFT(nSamples, nShorts, pRealIn, pRealOut, pImageOut);
// 对获得的FFT变换结果进行自乘，计算功率谱
for (unsigned int i = 0; i < nSamples; i++)
{
pDataOut[i] = sqrt(pRealOut[i] * pRealOut[i] + pImageOut[i] * pImageOut[i]);
}
delete[] pRealOut;
delete[] pImageOut;
}
}
//
// 获取输入信号的对数功率谱
void CWaveConvertor::ConvertToLogPowerSpectral(
unsigned int nSamples, // 样本数量
unsigned int nShorts, // 短时点数
const double* pRealIn, // 输入信号
double* pDataOut // 输出的对数功率谱
)
{
// 先转换为功率谱
CWaveConvertor::ConvertToPowerSpectral(nSamples, nShorts, pRealIn, pDataOut);
// 转化为对数功率谱
for (unsigned int i = 0; i < nSamples; i++)
{
pDataOut[i] = log(pDataOut[i]);
}
}
//
// 获取输入信号的倒谱
void CWaveConvertor::ConvertToCepStrum(
unsigned int nSamples, // 样本数量
unsigned int nShorts, // 短时点数
const double* pRealIn, // 输入信号
double* pDataOut // 输出的倒谱
)
{
// 保存对数功率谱
double* pTemp = new double[nSamples];
// 保存虚部数据
double* pImageOut = new double[nSamples];
// 获取对数功率谱
CWaveConvertor::ConvertToLogPowerSpectral(nSamples, nShorts, pRealIn, pTemp);
// 再次进行傅立叶逆变换，得到倒谱
CWaveConvertor::ConvertToRFT(nSamples, nShorts, pTemp, pDataOut, pImageOut);
delete[] pTemp;
delete[] pImageOut;
}
//
// 对输入的8位信号，转化为double类型序列
void CWaveConvertor::ConvertToDoubleMono(
const byte* pDataIn, // 输入样本序列
unsigned int nCount, // 样本数量
double* pDataOut // 输出double类型序列
)
{
int nAdjust = 127;
if (pDataIn != NULL)
{
// 填零
memset(pDataOut, 0, sizeof(double) * nCount);
// 8位信号减127
for (unsigned int i = 0; i < nCount; i++)
{
pDataOut[i] = pDataIn[i] - nAdjust;
}
}
}
//
// 对输入的16位信号，转化为double类型序列
void CWaveConvertor::ConvertToDoubleMono(
const int* pDataIn, // 输入样本序列
unsigned int nCount, // 样本数量
double* pDataOut // 输出double类型序列
)
{
char* pDataTemp = NULL;
unsigned int nOffset = 0;
if (pDataIn != NULL)
{
pDataTemp = (char*) pDataIn;
// 填零
nCount /= 2;
memset(pDataOut, 0, sizeof(double) * nCount);
for (unsigned int i = 0; i < nCount; i++)
{
nOffset = i * 2;
// 对高低位倒置的16位信号进行计算
pDataOut[i] = pDataTemp[nOffset] + pDataTemp[nOffset + 1] * 0x100;
}
}
}
//
// 对输入的8位信号，转化为双声道double序列
void CWaveConvertor::ConvertToDoubleStereo(
const byte* pDataIn, // 输入样本序列
unsigned int nCount, // 样本数量
double* pDataOutLeft, // 输出左声道double类型序列
double* pDataOutRight // 输出右声道double类型序列
)
{
double* pDataTemp = NULL;
if (pDataIn != NULL)
{
pDataTemp = new double[nCount];
// 信号类型先转化为double类型
CWaveConvertor::ConvertToDoubleMono(pDataIn, nCount, pDataTemp);
//提取左右声道
for (unsigned int i = 0; i < nCount; i += 2)
{
pDataOutLeft[i / 2] = pDataTemp[i];
pDataOutRight[i / 2] = pDataTemp[i + 1];
}
delete[] pDataTemp;
}
}
// 对输入的16位信号，转化为双声道double序列
void CWaveConvertor::ConvertToDoubleStereo(
const int* pDataIn, // 输入样本序列
unsigned int nCount, // 样本数量
double* pDataOutLeft, // 输出左声道double类型序列
double* pDataOutRight // 输出右声道double类型序列
)
{
double* pDataTemp = NULL;
if (pDataIn != NULL)
{
pDataTemp = new double[nCount];
// 信号类型先转化为double类型
CWaveConvertor::ConvertToDoubleMono(pDataIn, nCount, pDataTemp);
//提取左右声道
for (unsigned int i = 0; i < nCount; i += 2)
{
pDataOutLeft[i / 2] = pDataTemp[i];
pDataOutRight[i / 2] = pDataTemp[i + 1];
}
delete[] pDataTemp;
}
}
//
// 对输入信号进行抽样
void CWaveConvertor::ConvertToSample(
const double* pDataIn, // 输入样本序列
unsigned int nCount, // 样本数量
double* pDataOut, // 输出右声道double类型序列
unsigned int nWinSize // 抽样窗口宽度
)
{
double* pDataTemp = NULL;
unsigned int nOffSet = 0;
if (pDataIn != NULL)
{
pDataTemp = new double[nCount / nWinSize];
// 进行抽样
for (unsigned int i = 0; i < nCount; i += nWinSize)
{
pDataOut[nOffSet] = pDataIn[i];
nOffSet ++;
}
delete[] pDataTemp;
}
}
//
// 对输入信号进行截取
void CWaveConvertor::GetData(
const double* pDataIn, // 输入样本序列
unsigned int nCount, // 输入样本数量
unsigned int nStart, // 起始下标
unsigned int nLen, // 截取长度
double* pDataOut // 截取后的序列
)
{
if (pDataIn != NULL && nCount >= (nStart + nLen))
{
memcpy(pDataOut, pDataIn + nStart, sizeof(double) * nLen);
}
}
