WaveConvertor.cpp
上传用户:goak128
上传日期:2013-07-17
资源大小:155k
文件大小:15k
源码类别:

控制台编程

开发平台:

C/C++

  1. #include "StdAfx.h"
  2. #include ".waveconvertor.h"
  3. //////////////////////////////////////////////////////////////////////////
  4. // class CWaveConvertor
  5. //
  6. // 功能: 实现裸音频数据的FFT变换
  7. // 创建人: 陈文凯 (chwkai@gmail.com)
  8. // 创建日期:2005年5月19日
  9. // 修改人:
  10. // 修改日期:
  11. // 版本
  12. /* Pjotr '87.  */
  13. // As best as I can determine this is in the public domain. No license was found
  14. // withit what so ever. The author is listed as:
  15. // Peter Valkenburg (valke@cs.vu.nl).
  16. //
  17. // If there is some disagreement then contact Bruce Forsberg (forsberg@tns.net)
  18. #include <math.h>
  19. #define pi   3.1415926535897932384626434
  20. #define c_re(c) ((c).re)
  21. #define c_im(c) ((c).im)
  22. /* C_add_mul adds product of c1 and c2 to c.  */
  23. #define c_add_mul(c, c1, c2) { COMPLEX C1, C2; C1 = (c1); C2 = (c2); 
  24. c_re (c) += C1.re * C2.re - C1.im * C2.im; 
  25. c_im (c) += C1.re * C2.im + C1.im * C2.re; }
  26. /* C_conj substitutes c by its complex conjugate. */
  27. #define c_conj(c) { c_im (c) = -c_im (c); }
  28. /* C_realdiv divides complex c by real.  */
  29. #define c_realdiv(c, real) { c_re (c) /= (real); c_im (c) /= (real); }
  30. /*
  31. * W gives the (already computed) Wn ^ k (= e ^ (2pi * i * k / n)).
  32. * Notice that the powerseries of Wn has period Nfactors.
  33. */
  34. #define W(n, k) (W_factors [((k) * (Nfactors / (n))) % Nfactors])
  35. /*! brief Constructor.
  36. */
  37. aflibFFT::aflibFFT()
  38. {
  39. Nfactors = 0;
  40. W_factors = NULL;
  41. }
  42. /*! brief Destructor.
  43. */
  44. aflibFFT::~aflibFFT()
  45. {
  46. if (W_factors != NULL)
  47. delete W_factors;
  48. }
  49. /*! brief Performs a forward or reverse FFT.
  50. This is the main API is this class. It will perform either a forward or
  51. inverse FFT depending how InverseTransform is set. If set to FALSE then
  52. forward FFT will be performed, TRUE and a inverse FFT will be performed.
  53. The number of samlpes (NumSamples) must be a power of 2. The ImagIn
  54. pointer can be NULL if there are no imaginary values. The user is
  55. responsable for passing in pointers for RealOut and ImagOut containing
  56. arrays of the proper size.
  57. */
  58. void
  59. aflibFFT::fft_double (
  60.   unsigned  NumSamples,
  61.   int       InverseTransform,
  62.   const double   *RealIn,
  63.   const double   *ImagIn,
  64.   double   *RealOut,
  65.   double   *ImagOut )
  66. {
  67. COMPLEX  in[1024];
  68. COMPLEX  out[1024];
  69. COMPLEX  * in_local = NULL;
  70. COMPLEX  * out_local = NULL;
  71. register COMPLEX  * in_ptr;
  72. register COMPLEX  * out_ptr;
  73. register unsigned int      i;
  74. // IF 1024 samples or less use local buffer else allocate memory
  75. if (NumSamples > 1024)
  76. {
  77. in_local = new COMPLEX[NumSamples];
  78. out_local = new COMPLEX[NumSamples];
  79. in_ptr = in_local;
  80. out_ptr = out_local;
  81. }
  82. else
  83. {
  84. in_ptr = in;
  85. out_ptr = out;
  86. }
  87. // Fill real and imaginary array
  88. for (i = 0; i < NumSamples; i++)
  89. {
  90. c_re(in_ptr[i]) = RealIn[i];
  91. if (ImagIn == NULL)
  92. c_im(in_ptr[i]) = 0.0;
  93. else
  94. c_im(in_ptr[i]) = ImagIn[i];
  95. }
  96. // Perform transform
  97. if (InverseTransform == TRUE)
  98. {
  99. rft(in_ptr, NumSamples, out_ptr);
  100. }
  101. else
  102. {
  103. fft(in_ptr, NumSamples, out_ptr);
  104. }
  105. // Fill real and imaginary array
  106. for (i = 0; i < NumSamples; i++)
  107. {
  108. RealOut[i] = c_re(out_ptr[i]);
  109. ImagOut[i] = c_im(out_ptr[i]);
  110. }
  111. // Free memory if local arrays were not used
  112. if (in_local != NULL)
  113. delete [] in_local;
  114. if (out_local != NULL)
  115. delete [] out_local;
  116. }
  117. /*
  118. * Forward Fast Fourier Transform on the n samples of complex array in.
  119. * The result is placed in out.  The number of samples, n, is arbitrary.
  120. * The W-factors are calculated in advance.
  121. */
  122. int
  123. aflibFFT::fft (
  124.    COMPLEX *in,
  125.    unsigned  n,
  126.    COMPLEX *out)
  127. {
  128. unsigned i;
  129. for (i = 0; i < n; i++)
  130. c_conj (in [i]);
  131. if (W_init (n) == -1)
  132. return -1;
  133. Fourier (in, n, out);
  134. for (i = 0; i < n; i++) {
  135. c_conj (out [i]);
  136. c_realdiv (out [i], n);
  137. }
  138. return 0;
  139. }
  140. /*
  141. * Reverse Fast Fourier Transform on the n complex samples of array in.
  142. * The result is placed in out.  The number of samples, n, is arbitrary.
  143. * The W-factors are calculated in advance.
  144. */
  145. int
  146. aflibFFT::rft (
  147.    COMPLEX *in,
  148.    unsigned  n,
  149.    COMPLEX *out)
  150. {
  151. if (W_init (n) == -1)
  152. return -1;
  153. Fourier (in, n, out);
  154. return 0;
  155. }
  156. /*
  157. * Recursive (reverse) complex fast Fourier transform on the n
  158. * complex samples of array in, with the Cooley-Tukey method.
  159. * The result is placed in out.  The number of samples, n, is arbitrary.
  160. * The algorithm costs O (n * (r1 + .. + rk)), where k is the number
  161. * of factors in the prime-decomposition of n (also the maximum
  162. * depth of the recursion), and ri is the i-th primefactor.
  163. */
  164. void
  165. aflibFFT::Fourier (
  166.    COMPLEX *in,
  167.    unsigned  n,
  168.    COMPLEX *out)
  169. {
  170. unsigned r;
  171. if ((r = radix (n)) < n)
  172. split (in, r, n / r, out);
  173. join (in, n / r, n, out);
  174. }
  175. /*
  176. * Give smallest possible radix for n samples.
  177. * Determines (in a rude way) the smallest primefactor of n.
  178. */
  179. unsigned
  180. aflibFFT::radix (unsigned n)
  181. {
  182. unsigned r;
  183. if (n < 2)
  184. return 1;
  185. for (r = 2; r < n; r++)
  186. if (n % r == 0)
  187. break;
  188. return r;
  189. }
  190. /*
  191. * Split array in of r * m samples in r parts of each m samples,
  192. * such that in [i] goes to out [(i % r) * m + (i / r)].
  193. * Then call for each part of out Fourier, so the r recursively
  194. * transformed parts will go back to in.
  195. */
  196. void
  197. aflibFFT::split (
  198.  register COMPLEX *in,
  199.  register unsigned r,
  200.  register unsigned m,
  201.  register COMPLEX *out)
  202. {
  203. register unsigned k, s, i, j;
  204. for (k = 0, j = 0; k < r; k++)
  205. for (s = 0, i = k; s < m; s++, i += r, j++)
  206. out [j] = in [i];
  207. for (k = 0; k < r; k++, out += m, in += m)
  208. Fourier (out, m, in);
  209. }
  210. /*
  211. * Sum the n / m parts of each m samples of in to n samples in out.
  212. *     r - 1
  213. * Out [j] becomes  sum  in [j % m] * W (j * k).  Here in is the k-th
  214. *     k = 0   k        n  k
  215. * part of in (indices k * m ... (k + 1) * m - 1), and r is the radix.
  216. * For k = 0, a complex multiplication with W (0) is avoided.
  217. */
  218. void
  219. aflibFFT::join (
  220. register COMPLEX *in,
  221. register unsigned m,
  222. register unsigned n,
  223. register COMPLEX *out)
  224. {
  225. register unsigned i, j, jk, s;
  226. for (s = 0; s < m; s++)
  227. for (j = s; j < n; j += m) {
  228. out [j] = in [s];
  229. for (i = s + m, jk = j; i < n; i += m, jk += j)
  230. c_add_mul (out [j], in [i], W (n, jk));
  231. }
  232. }
  233. int
  234. aflibFFT::W_init(unsigned n)
  235. {
  236. unsigned k;
  237. if (n == Nfactors)
  238. return 0;
  239. if (Nfactors != 0 && W_factors != 0)
  240. delete [] W_factors;
  241. if ((Nfactors = n) == 0)
  242. return 0;
  243. if ((W_factors = new COMPLEX[n]) == NULL)
  244. return -1;
  245. for (k = 0; k < n; k++) {
  246. c_re (W_factors [k]) = cos (2 * pi * k / n);
  247. c_im (W_factors [k]) = sin (2 * pi * k / n);
  248. }
  249. return 0;
  250. }
  251. //////////////////////////////////////////////////////////////////////////
  252. // 构造函数
  253. CWaveConvertor::CWaveConvertor(void)
  254. {
  255. }
  256. //////////////////////////////////////////////////////////////////////////
  257. // 析构函数
  258. CWaveConvertor::~CWaveConvertor(void)
  259. {
  260. }
  261. //////////////////////////////////////////////////////////////////////////
  262. // 对输入数据进行短时快速FFT变换,输入数据为已经加窗的数据
  263. void CWaveConvertor::ConverToFFT( 
  264.  unsigned int nSamples, // 样本数量
  265.  unsigned int nShorts, // 短时点数
  266.  const double* pRealIn, // 输入的实部,不可为空
  267.  double* pRealOut, // 输出的实部
  268.  double* pImageOut // 输出的虚部
  269.  )
  270. {
  271. aflibFFT libFFT;
  272. // 计算所需变换次数
  273. unsigned int nCount = nSamples / nShorts;
  274. // 数组偏移
  275. unsigned int nOffSet = 0;
  276. // 初始化输出数据
  277. memset(pRealOut, 0, sizeof(double) * nSamples);
  278. memset(pImageOut, 0, sizeof(double) * nSamples);
  279. for (unsigned int i = 0; i < nCount; i++)
  280. {
  281. nOffSet = i * nShorts;
  282. // 进行快速FFT
  283. libFFT.fft_double(
  284. nShorts, 
  285. FALSE, 
  286. &pRealIn[nOffSet], 
  287. NULL, 
  288. &pRealOut[nOffSet], 
  289. &pImageOut[nOffSet]
  290. );
  291. }
  292. }
  293. //////////////////////////////////////////////////////////////////////////
  294. // 对输入数据进行Reverse Fourier转化,输入数据为已经加窗的数据
  295. void CWaveConvertor::ConvertToRFT( 
  296.   unsigned int nSamples, // 样本数量
  297.   unsigned int nShorts, // 短时点数
  298.   const double* pRealIn, // 输入的实部,不可为空
  299.   double* pRealOut, // 输出的实部
  300.   double* pImageOut // 输出的虚部
  301.   )
  302. {
  303. aflibFFT libFFT;
  304. // 计算所需变换次数
  305. unsigned int nCount = nSamples / nShorts;
  306. // 数组偏移
  307. unsigned int nOffSet = 0;
  308. // 初始化输出数据
  309. memset(pRealOut, 0, sizeof(double) * nSamples);
  310. memset(pImageOut, 0, sizeof(double) * nSamples);
  311. for (unsigned int i = 0; i < nCount; i++)
  312. {
  313. nOffSet = i * nShorts;
  314. // 进行快速FFT
  315. libFFT.fft_double(
  316. nShorts, 
  317. TRUE, 
  318. &pRealIn[nOffSet], 
  319. NULL, 
  320. &pRealOut[nOffSet], 
  321. &pImageOut[nOffSet]
  322. );
  323. }
  324. }
  325. //////////////////////////////////////////////////////////////////////////
  326. // 获取输入信号的功率谱
  327. void CWaveConvertor::ConvertToPowerSpectral( 
  328. unsigned int nSamples, // 样本数量
  329. unsigned int nShorts, // 短时点数
  330. const double* pRealIn, // 输入信号
  331. double* pDataOut // 输出的功率谱
  332. )
  333. {
  334. double* pRealOut = NULL;
  335. double* pImageOut = NULL;
  336. // 输入数据是否有效
  337. if (nSamples > 0 && pRealIn != NULL)
  338. {
  339. // 为FFT输出数据分配空间
  340. pRealOut = new double[nSamples];
  341. pImageOut = new double[nSamples];
  342. // 进行FFT变换
  343. CWaveConvertor::ConverToFFT(nSamples, nShorts, pRealIn, pRealOut, pImageOut);
  344. // 对获得的FFT变换结果进行自乘,计算功率谱
  345. for (unsigned int i = 0; i < nSamples; i++)
  346. {
  347. pDataOut[i] = sqrt(pRealOut[i] * pRealOut[i] + pImageOut[i] * pImageOut[i]);
  348. }
  349. delete[] pRealOut;
  350. delete[] pImageOut;
  351. }
  352. }
  353. //////////////////////////////////////////////////////////////////////////
  354. // 获取输入信号的对数功率谱
  355. void CWaveConvertor::ConvertToLogPowerSpectral( 
  356. unsigned int nSamples, // 样本数量
  357. unsigned int nShorts, // 短时点数
  358. const double* pRealIn, // 输入信号
  359. double* pDataOut // 输出的对数功率谱
  360. )
  361. {
  362. // 先转换为功率谱
  363. CWaveConvertor::ConvertToPowerSpectral(nSamples, nShorts, pRealIn, pDataOut);
  364. // 转化为对数功率谱
  365. for (unsigned int i = 0; i < nSamples; i++)
  366. {
  367. pDataOut[i] = log(pDataOut[i]);
  368. }
  369. }
  370. //////////////////////////////////////////////////////////////////////////
  371. // 获取输入信号的倒谱
  372. void CWaveConvertor::ConvertToCepStrum( 
  373.    unsigned int nSamples, // 样本数量
  374.    unsigned int nShorts, // 短时点数
  375.    const double* pRealIn, // 输入信号
  376.    double* pDataOut // 输出的倒谱
  377.    )
  378. {
  379. // 保存对数功率谱
  380. double* pTemp = new double[nSamples];
  381. // 保存虚部数据
  382. double* pImageOut = new double[nSamples];
  383. // 获取对数功率谱
  384. CWaveConvertor::ConvertToLogPowerSpectral(nSamples, nShorts, pRealIn, pTemp);
  385. // 再次进行傅立叶逆变换,得到倒谱
  386. CWaveConvertor::ConvertToRFT(nSamples, nShorts, pTemp, pDataOut, pImageOut);
  387. delete[] pTemp;
  388. delete[] pImageOut;
  389. }
  390. //////////////////////////////////////////////////////////////////////////
  391. // 对输入的8位信号,转化为double类型序列
  392. void CWaveConvertor::ConvertToDoubleMono( 
  393. const byte* pDataIn, // 输入样本序列
  394. unsigned int nCount, // 样本数量
  395. double* pDataOut // 输出double类型序列
  396. )
  397. {
  398. int nAdjust = 127;
  399. if (pDataIn != NULL)
  400. {
  401. // 填零
  402. memset(pDataOut, 0, sizeof(double) * nCount);
  403. // 8位信号减127
  404. for (unsigned int i = 0; i < nCount; i++)
  405. {
  406. pDataOut[i] = pDataIn[i] - nAdjust;
  407. }
  408. }
  409. }
  410. //////////////////////////////////////////////////////////////////////////
  411. // 对输入的16位信号,转化为double类型序列
  412. void CWaveConvertor::ConvertToDoubleMono( 
  413. const int* pDataIn, // 输入样本序列
  414. unsigned int nCount, // 样本数量
  415. double* pDataOut // 输出double类型序列
  416. )
  417. {
  418. char* pDataTemp = NULL;
  419. unsigned int nOffset = 0;
  420. if (pDataIn != NULL)
  421. {
  422. pDataTemp = (char*) pDataIn;
  423. // 填零
  424. nCount /= 2;
  425. memset(pDataOut, 0, sizeof(double) * nCount);
  426. for (unsigned int i = 0; i < nCount; i++)
  427. {
  428. nOffset = i * 2;
  429. // 对高低位倒置的16位信号进行计算
  430. pDataOut[i] = pDataTemp[nOffset] + pDataTemp[nOffset + 1] * 0x100;
  431. }
  432. }
  433. }
  434. //////////////////////////////////////////////////////////////////////////
  435. // 对输入的8位信号,转化为双声道double序列
  436. void CWaveConvertor::ConvertToDoubleStereo(
  437. const byte* pDataIn, // 输入样本序列
  438. unsigned int nCount, // 样本数量
  439. double* pDataOutLeft, // 输出左声道double类型序列
  440. double* pDataOutRight // 输出右声道double类型序列
  441. )
  442. {
  443. double* pDataTemp = NULL;
  444. if (pDataIn != NULL)
  445. {
  446. pDataTemp = new double[nCount];
  447. // 信号类型先转化为double类型
  448. CWaveConvertor::ConvertToDoubleMono(pDataIn, nCount, pDataTemp);
  449. //提取左右声道
  450. for (unsigned int i = 0; i < nCount; i += 2)
  451. {
  452. pDataOutLeft[i / 2] = pDataTemp[i];
  453. pDataOutRight[i / 2] = pDataTemp[i + 1];
  454. }
  455. delete[] pDataTemp;
  456. }
  457. }
  458. // 对输入的16位信号,转化为双声道double序列
  459. void CWaveConvertor::ConvertToDoubleStereo(
  460. const int* pDataIn, // 输入样本序列
  461. unsigned int nCount, // 样本数量
  462. double* pDataOutLeft, // 输出左声道double类型序列
  463. double* pDataOutRight // 输出右声道double类型序列
  464. )
  465. {
  466. double* pDataTemp = NULL;
  467. if (pDataIn != NULL)
  468. {
  469. pDataTemp = new double[nCount];
  470. // 信号类型先转化为double类型
  471. CWaveConvertor::ConvertToDoubleMono(pDataIn, nCount, pDataTemp);
  472. //提取左右声道
  473. for (unsigned int i = 0; i < nCount; i += 2)
  474. {
  475. pDataOutLeft[i / 2] = pDataTemp[i];
  476. pDataOutRight[i / 2] = pDataTemp[i + 1];
  477. }
  478. delete[] pDataTemp;
  479. }
  480. }
  481. //////////////////////////////////////////////////////////////////////////
  482. // 对输入信号进行抽样
  483. void CWaveConvertor::ConvertToSample(
  484. const double* pDataIn, // 输入样本序列
  485. unsigned int nCount, // 样本数量
  486. double* pDataOut, // 输出右声道double类型序列
  487. unsigned int nWinSize // 抽样窗口宽度
  488. )
  489. {
  490. double* pDataTemp = NULL;
  491. unsigned int nOffSet = 0;
  492. if (pDataIn != NULL)
  493. {
  494. pDataTemp = new double[nCount / nWinSize];
  495. // 进行抽样
  496. for (unsigned int i = 0; i < nCount; i += nWinSize)
  497. {
  498. pDataOut[nOffSet] = pDataIn[i];
  499. nOffSet ++;
  500. }
  501. delete[] pDataTemp;
  502. }
  503. }
  504. //////////////////////////////////////////////////////////////////////////
  505. // 对输入信号进行截取
  506. void CWaveConvertor::GetData(
  507.  const double* pDataIn, // 输入样本序列
  508.  unsigned int nCount, // 输入样本数量
  509.  unsigned int nStart, // 起始下标
  510.  unsigned int nLen, // 截取长度
  511.  double* pDataOut // 截取后的序列
  512.  )
  513. {
  514. if (pDataIn != NULL && nCount >= (nStart + nLen))
  515. {
  516. memcpy(pDataOut, pDataIn + nStart, sizeof(double) * nLen);
  517. }
  518. }