开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 流媒体/Mpeg4/MP4 > 查看源码
filtbank.c
资源名称：NETVIDEO.rar　[点击查看]
上传用户：sun1608
上传日期：2007-02-02
资源大小：6116k
文件大小：18k
源码类别：

流媒体/Mpeg4/MP4

开发平台：

Visual C++

filtbank.c：源码内容
  1. /*
  2.  * FAAC - Freeware Advanced Audio Coder
  3.  * Copyright (C) 2001 Menno Bakker
  4.  *
  5.  * This library is free software; you can redistribute it and/or
  6.  * modify it under the terms of the GNU Lesser General Public
  7.  * License as published by the Free Software Foundation; either
  8.  * version 2.1 of the License, or (at your option) any later version.
  9.  *
  10.  * This library is distributed in the hope that it will be useful,
  11.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * Lesser General Public License for more details.
  15.  * You should have received a copy of the GNU Lesser General Public
  16.  * License along with this library; if not, write to the Free Software
  17.  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18.  *
  19.  * $Id: filtbank.c,v 1.3 2001/06/04 23:02:24 wmay Exp $
  20.  */
  22. /*
  23.  * CHANGES:
  24.  *  2001/01/17: menno: Added frequency cut off filter.
  25.  *
  26.  */
  28. #include <math.h>
  29. #include <stdio.h>
  30. #include <stdlib.h>
  32. #include "coder.h"
  33. #include "filtbank.h"
  34. #include "frame.h"
  35. #include "fft.h"
  36. #include "util.h"
  38. #define  TWOPI       2*M_PI
  41. void FilterBankInit(faacEncHandle hEncoder)
  42. {
  43. unsigned int i, channel;
  45. for (channel = 0; channel < hEncoder->numChannels; channel++) {
  46. hEncoder->freqBuff[channel] = (double*)AllocMemory(2*FRAME_LEN*sizeof(double));
  47. hEncoder->overlapBuff[channel] = (double*)AllocMemory(FRAME_LEN*sizeof(double));
  48. SetMemory(hEncoder->overlapBuff[channel], 0, FRAME_LEN*sizeof(double));
  49. }
  51. hEncoder->sin_window_long = (double*)AllocMemory(BLOCK_LEN_LONG*sizeof(double));
  52. hEncoder->sin_window_short = (double*)AllocMemory(BLOCK_LEN_SHORT*sizeof(double));
  53. hEncoder->kbd_window_long = (double*)AllocMemory(BLOCK_LEN_LONG*sizeof(double));
  54. hEncoder->kbd_window_short = (double*)AllocMemory(BLOCK_LEN_SHORT*sizeof(double));
  56. for( i=0; i<BLOCK_LEN_LONG; i++ )
  57. hEncoder->sin_window_long[i] = sin((M_PI/(2*BLOCK_LEN_LONG)) * (i + 0.5));
  58. for( i=0; i<BLOCK_LEN_SHORT; i++ )
  59. hEncoder->sin_window_short[i] = sin((M_PI/(2*BLOCK_LEN_SHORT)) * (i + 0.5));
  61. CalculateKBDWindow(hEncoder->kbd_window_long, 4, BLOCK_LEN_LONG*2);
  62. CalculateKBDWindow(hEncoder->kbd_window_short, 6, BLOCK_LEN_SHORT*2);
  63. }
  65. void FilterBankEnd(faacEncHandle hEncoder)
  66. {
  67. unsigned int channel;
  69. for (channel = 0; channel < hEncoder->numChannels; channel++) {
  70. if (hEncoder->freqBuff[channel]) FreeMemory(hEncoder->freqBuff[channel]);
  71. if (hEncoder->overlapBuff[channel]) FreeMemory(hEncoder->overlapBuff[channel]);
  72. }
  74. if (hEncoder->sin_window_long) FreeMemory(hEncoder->sin_window_long);
  75. if (hEncoder->sin_window_short) FreeMemory(hEncoder->sin_window_short);
  76. if (hEncoder->kbd_window_long) FreeMemory(hEncoder->kbd_window_long);
  77. if (hEncoder->kbd_window_short) FreeMemory(hEncoder->kbd_window_short);
  78. }
  80. void FilterBank(faacEncHandle hEncoder,
  81. CoderInfo *coderInfo,
  82. double *p_in_data,
  83. double *p_out_mdct,
  84. double *p_overlap,
  85. int overlap_select)
  86. {
  87. double *p_o_buf, *first_window, *second_window;
  88. double *transf_buf;
  89. int k, i;
  90. int block_type = coderInfo->block_type;
  92. transf_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
  94. /* create / shift old values */
  95. /* We use p_overlap here as buffer holding the last frame time signal*/
  96. if(overlap_select != MNON_OVERLAPPED) {
  97. memcpy(transf_buf, p_overlap, FRAME_LEN*sizeof(double));
  98. memcpy(transf_buf+BLOCK_LEN_LONG, p_in_data, FRAME_LEN*sizeof(double));
  99. memcpy(p_overlap, p_in_data, FRAME_LEN*sizeof(double));
  100. } else {
  101. memcpy(transf_buf, p_in_data, 2*FRAME_LEN*sizeof(double));
  102. }
  104. /*  Window shape processing */
  105. if(overlap_select != MNON_OVERLAPPED) {
  106. switch (coderInfo->prev_window_shape) {
  107. case SINE_WINDOW:
  108. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
  109. first_window = hEncoder->sin_window_long;
  110. else
  111. first_window = hEncoder->sin_window_short;
  112. break;
  113. case KBD_WINDOW:
  114. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
  115. first_window = hEncoder->kbd_window_long;
  116. else
  117. first_window = hEncoder->kbd_window_short;
  118. break;
  119. }
  121. switch (coderInfo->window_shape){
  122. case SINE_WINDOW:
  123. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
  124. second_window = hEncoder->sin_window_long;
  125. else
  126. second_window = hEncoder->sin_window_short;
  127. break;
  128. case KBD_WINDOW:
  129. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
  130. second_window = hEncoder->kbd_window_long;
  131. else
  132. second_window = hEncoder->kbd_window_short;
  133. break;
  134. }
  135. } else {
  136. /* Always long block and sine window for LTP */
  137. first_window = hEncoder->sin_window_long;
  138. second_window = hEncoder->sin_window_long;
  139. }
  141. /* Set ptr to transf-Buffer */
  142. p_o_buf = transf_buf;
  144. /* Separate action for each Block Type */
  145. switch (block_type) {
  146.     case ONLY_LONG_WINDOW :
  147. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++){
  148. p_out_mdct[i] = p_o_buf[i] * first_window[i];
  149. p_out_mdct[i+BLOCK_LEN_LONG] = p_o_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
  150. }
  151. MDCT( p_out_mdct, 2*BLOCK_LEN_LONG );
  152. break;
  154.     case LONG_SHORT_WINDOW :
  155. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  156. p_out_mdct[i] = p_o_buf[i] * first_window[i];
  157. memcpy(p_out_mdct+BLOCK_LEN_LONG,p_o_buf+BLOCK_LEN_LONG,NFLAT_LS*sizeof(double));
  158. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
  159. p_out_mdct[i+BLOCK_LEN_LONG+NFLAT_LS] = p_o_buf[i+BLOCK_LEN_LONG+NFLAT_LS] * second_window[BLOCK_LEN_SHORT-i-1];
  160. SetMemory(p_out_mdct+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
  161. MDCT( p_out_mdct, 2*BLOCK_LEN_LONG );
  162. break;
  164.     case SHORT_LONG_WINDOW :
  165. SetMemory(p_out_mdct,0,NFLAT_LS*sizeof(double));
  166. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
  167. p_out_mdct[i+NFLAT_LS] = p_o_buf[i+NFLAT_LS] * first_window[i];
  168. memcpy(p_out_mdct+NFLAT_LS+BLOCK_LEN_SHORT,p_o_buf+NFLAT_LS+BLOCK_LEN_SHORT,NFLAT_LS*sizeof(double));
  169. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  170. p_out_mdct[i+BLOCK_LEN_LONG] = p_o_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
  171. MDCT( p_out_mdct, 2*BLOCK_LEN_LONG );
  172. break;
  174.     case ONLY_SHORT_WINDOW :
  175. p_o_buf += NFLAT_LS;
  176. for ( k=0; k < MAX_SHORT_WINDOWS; k++ ) {
  177. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++ ){
  178. p_out_mdct[i] = p_o_buf[i] * first_window[i];
  179. p_out_mdct[i+BLOCK_LEN_SHORT] = p_o_buf[i+BLOCK_LEN_SHORT] * second_window[BLOCK_LEN_SHORT-i-1];
  180. }
  181. MDCT( p_out_mdct, 2*BLOCK_LEN_SHORT );
  182. p_out_mdct += BLOCK_LEN_SHORT;
  183. p_o_buf += BLOCK_LEN_SHORT;
  184. first_window = second_window;
  185. }
  186. break;
  187. }
  189. if (transf_buf) FreeMemory(transf_buf);
  190. }
  192. void IFilterBank(faacEncHandle hEncoder,
  193.  CoderInfo *coderInfo,
  194.  double *p_in_data,
  195.  double *p_out_data,
  196.  double *p_overlap,
  197.  int overlap_select)
  198. {
  199. double *o_buf, *transf_buf, *overlap_buf;
  200. double *first_window, *second_window;
  202. double  *fp;
  203. int k, i;
  204. int block_type = coderInfo->block_type;
  206. transf_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
  207. overlap_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
  209. /*  Window shape processing */
  210. if (overlap_select != MNON_OVERLAPPED) {
  211. // switch (coderInfo->prev_window_shape){
  212. // case SINE_WINDOW:
  213. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
  214. first_window = hEncoder->sin_window_long;
  215. else
  216. first_window = hEncoder->sin_window_short;
  217. // break;
  218. // case KBD_WINDOW:
  219. // if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
  220. // first_window = hEncoder->kbd_window_long;
  221. // else
  222. // first_window = hEncoder->kbd_window_short;
  223. // break;
  224. // }
  226. // switch (coderInfo->window_shape){
  227. // case SINE_WINDOW:
  228. if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
  229. second_window = hEncoder->sin_window_long;
  230. else
  231. second_window = hEncoder->sin_window_short;
  232. // break;
  233. // case KBD_WINDOW:
  234. // if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
  235. // second_window = hEncoder->kbd_window_long;
  236. // else
  237. // second_window = hEncoder->kbd_window_short;
  238. // break;
  239. // }
  240. } else {
  241. /* Always long block and sine window for LTP */
  242. first_window  = hEncoder->sin_window_long;
  243. second_window = hEncoder->sin_window_long;
  244. }
  246. /* Assemble overlap buffer */
  247. memcpy(overlap_buf,p_overlap,BLOCK_LEN_LONG*sizeof(double));
  248. o_buf = overlap_buf;
  250. /* Separate action for each Block Type */
  251. switch( block_type ) {
  252.     case ONLY_LONG_WINDOW :
  253. memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
  254. IMDCT( transf_buf, 2*BLOCK_LEN_LONG );
  255. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  256. transf_buf[i] *= first_window[i];
  257. if (overlap_select != MNON_OVERLAPPED) {
  258. for ( i = 0 ; i < BLOCK_LEN_LONG; i++ ){
  259. o_buf[i] += transf_buf[i];
  260. o_buf[i+BLOCK_LEN_LONG] = transf_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
  261. }
  262. } else { /* overlap_select == NON_OVERLAPPED */
  263. for ( i = 0 ; i < BLOCK_LEN_LONG; i++ )
  264. transf_buf[i+BLOCK_LEN_LONG] *= second_window[BLOCK_LEN_LONG-i-1];
  265.         }
  266. break;
  268.     case LONG_SHORT_WINDOW :
  269. memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
  270. IMDCT( transf_buf, 2*BLOCK_LEN_LONG );
  271. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  272. transf_buf[i] *= first_window[i];
  273.         if (overlap_select != MNON_OVERLAPPED) {
  274. for ( i = 0 ; i < BLOCK_LEN_LONG; i++ )
  275. o_buf[i] += transf_buf[i];
  276. memcpy(o_buf+BLOCK_LEN_LONG,transf_buf+BLOCK_LEN_LONG,NFLAT_LS*sizeof(double));
  277. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
  278. o_buf[i+BLOCK_LEN_LONG+NFLAT_LS] = transf_buf[i+BLOCK_LEN_LONG+NFLAT_LS] * second_window[BLOCK_LEN_SHORT-i-1];
  279. SetMemory(o_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
  280.         } else { /* overlap_select == NON_OVERLAPPED */
  281. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
  282. transf_buf[i+BLOCK_LEN_LONG+NFLAT_LS] *= second_window[BLOCK_LEN_SHORT-i-1];
  283. SetMemory(transf_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
  284. }
  285. break;
  287. case SHORT_LONG_WINDOW :
  288. memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
  289. IMDCT( transf_buf, 2*BLOCK_LEN_LONG );
  290. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
  291. transf_buf[i+NFLAT_LS] *= first_window[i];
  292. if (overlap_select != MNON_OVERLAPPED) {
  293. for ( i = 0 ; i < BLOCK_LEN_SHORT; i++ )
  294. o_buf[i+NFLAT_LS] += transf_buf[i+NFLAT_LS];
  295. memcpy(o_buf+BLOCK_LEN_SHORT+NFLAT_LS,transf_buf+BLOCK_LEN_SHORT+NFLAT_LS,NFLAT_LS*sizeof(double));
  296. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  297. o_buf[i+BLOCK_LEN_LONG] = transf_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
  298. } else { /* overlap_select == NON_OVERLAPPED */
  299. SetMemory(transf_buf,0,NFLAT_LS*sizeof(double));
  300. for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
  301. transf_buf[i+BLOCK_LEN_LONG] *= second_window[BLOCK_LEN_LONG-i-1];
  302. }
  303. break;
  305. case ONLY_SHORT_WINDOW :
  306. if (overlap_select != MNON_OVERLAPPED) {
  307. fp = o_buf + NFLAT_LS;
  308. } else { /* overlap_select == NON_OVERLAPPED */
  309. fp = transf_buf;
  310. }
  311. for ( k=0; k < MAX_SHORT_WINDOWS; k++ ) {
  312. memcpy(transf_buf,p_in_data,BLOCK_LEN_SHORT*sizeof(double));
  313. IMDCT( transf_buf, 2*BLOCK_LEN_SHORT );
  314. p_in_data += BLOCK_LEN_SHORT;
  315. if (overlap_select != MNON_OVERLAPPED) {
  316. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
  317. transf_buf[i] *= first_window[i];
  318. fp[i] += transf_buf[i];
  319. fp[i+BLOCK_LEN_SHORT] = transf_buf[i+BLOCK_LEN_SHORT] * second_window[BLOCK_LEN_SHORT-i-1];
  320. }
  321. fp += BLOCK_LEN_SHORT;
  322. } else { /* overlap_select == NON_OVERLAPPED */
  323. for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
  324. fp[i] *= first_window[i];
  325. fp[i+BLOCK_LEN_SHORT] *= second_window[BLOCK_LEN_SHORT-i-1];
  326. }
  327. fp += 2*BLOCK_LEN_SHORT;
  328. }
  329. first_window = second_window;
  330. }
  331. SetMemory(o_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
  332. break;
  333. }
  335. if (overlap_select != MNON_OVERLAPPED)
  336. memcpy(p_out_data,o_buf,BLOCK_LEN_LONG*sizeof(double));
  337. else  /* overlap_select == NON_OVERLAPPED */
  338. memcpy(p_out_data,transf_buf,2*BLOCK_LEN_LONG*sizeof(double));
  340. /* save unused output data */
  341. memcpy(p_overlap,o_buf+BLOCK_LEN_LONG,BLOCK_LEN_LONG*sizeof(double));
  343. if (overlap_buf) FreeMemory(overlap_buf);
  344. if (transf_buf) FreeMemory(transf_buf);
  345. }
  347. void specFilter(double *freqBuff,
  348. int sampleRate,
  349. int lowpassFreq,
  350. int specLen
  351. )
  352. {
  353. int lowpass,xlowpass;
  355. /* calculate the last line which is not zero */
  356. lowpass = (lowpassFreq * specLen) / (sampleRate>>1) + 1;
  357. xlowpass = (lowpass < specLen) ? lowpass : specLen ;
  359. SetMemory(freqBuff+xlowpass,0,(specLen-xlowpass)*sizeof(double));
  360. }
  362. static double Izero(double x)
  363. {
  364. const double IzeroEPSILON = 1E-41;  /* Max error acceptable in Izero */
  365. double sum, u, halfx, temp;
  366. int n;
  368. sum = u = n = 1;
  369. halfx = x/2.0;
  370. do {
  371. temp = halfx/(double)n;
  372. n += 1;
  373. temp *= temp;
  374. u *= temp;
  375. sum += u;
  376. } while (u >= IzeroEPSILON*sum);
  378. return(sum);
  379. }
  381. static void CalculateKBDWindow(double* win, double alpha, int length)
  382. {
  383. int i;
  384. double IBeta;
  385. double tmp;
  386. double sum = 0.0;
  388. alpha *= M_PI;
  389. IBeta = 1.0/Izero(alpha);
  391. /* calculate lower half of Kaiser Bessel window */
  392. for(i=0; i<(length>>1); i++) {
  393. tmp = 4.0*(double)i/(double)length - 1.0;
  394. win[i] = Izero(alpha*sqrt(1.0-tmp*tmp))*IBeta;
  395. sum += win[i];
  396. }
  398. sum = 1.0/sum;
  399. tmp = 0.0;
  401. /* calculate lower half of window */
  402. for(i=0; i<(length>>1); i++) {
  403. tmp += win[i];
  404. win[i] = sqrt(tmp*sum);
  405. }
  406. }
  408. static void MDCT(double *data, int N)
  409. {
  410. double *xi, *xr;
  411. double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
  412. double freq = TWOPI / N;
  413. double cosfreq8, sinfreq8;
  414. int i, n;
  416. xi = (double*)AllocMemory((N >> 2)*sizeof(double));
  417. xr = (double*)AllocMemory((N >> 2)*sizeof(double));
  419. /* prepare for recurrence relation in pre-twiddle */
  420. cfreq = cos (freq);
  421. sfreq = sin (freq);
  422. cosfreq8 = cos (freq * 0.125);
  423. sinfreq8 = sin (freq * 0.125);
  424. c = cosfreq8;
  425. s = sinfreq8;
  427. for (i = 0; i < (N >> 2); i++) {
  428. /* calculate real and imaginary parts of g(n) or G(p) */
  429. n = (N >> 1) - 1 - 2 * i;
  431. if (i < (N >> 3))
  432. tempr = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n]; /* use second form of e(n) for n = N / 2 - 1 - 2i */
  433. else
  434. tempr = data [(N >> 2) + n] - data [(N >> 2) - 1 - n]; /* use first form of e(n) for n = N / 2 - 1 - 2i */
  436. n = 2 * i;
  437. if (i < (N >> 3))
  438. tempi = data [(N >> 2) + n] - data [(N >> 2) - 1 - n]; /* use first form of e(n) for n=2i */
  439. else
  440. tempi = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n]; /* use second form of e(n) for n=2i*/
  442. /* calculate pre-twiddled FFT input */
  443. xr[i] = tempr * c + tempi * s;
  444. xi[i] = tempi * c - tempr * s;
  446. /* use recurrence to prepare cosine and sine for next value of i */
  447. cold = c;
  448. c = c * cfreq - s * sfreq;
  449. s = s * cfreq + cold * sfreq;
  450. }
  452. /* Perform in-place complex FFT of length N/4 */
  453. switch (N) {
  454.     case 256:
  455. srfft(xr, xi, 6);
  456. break;
  457.     case 2048:
  458. srfft(xr, xi, 9);
  459. }
  461. /* prepare for recurrence relations in post-twiddle */
  462. c = cosfreq8;
  463. s = sinfreq8;
  465. /* post-twiddle FFT output and then get output data */
  466. for (i = 0; i < (N >> 2); i++) {
  467. /* get post-twiddled FFT output  */
  468. tempr = 2. * (xr[i] * c + xi[i] * s);
  469. tempi = 2. * (xi[i] * c - xr[i] * s);
  471. /* fill in output values */
  472. data [2 * i] = -tempr;   /* first half even */
  473. data [(N >> 1) - 1 - 2 * i] = tempi;  /* first half odd */
  474. data [(N >> 1) + 2 * i] = -tempi;  /* second half even */
  475. data [N - 1 - 2 * i] = tempr;  /* second half odd */
  477. /* use recurrence to prepare cosine and sine for next value of i */
  478. cold = c;
  479. c = c * cfreq - s * sfreq;
  480. s = s * cfreq + cold * sfreq;
  481. }
  483. if (xr) FreeMemory(xr);
  484. if (xi) FreeMemory(xi);
  485. }
  487. static void IMDCT(double *data, int N)
  488. {
  489. double *xi, *xr;
  490. double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
  491. double freq = 2.0 * M_PI / N;
  492. double fac, cosfreq8, sinfreq8;
  493. int i;
  495. xi = (double*)AllocMemory((N >> 2)*sizeof(double));
  496. xr = (double*)AllocMemory((N >> 2)*sizeof(double));
  498. /* Choosing to allocate 2/N factor to Inverse Xform! */
  499. fac = 2. / N; /* remaining 2/N from 4/N IFFT factor */
  501. /* prepare for recurrence relation in pre-twiddle */
  502. cfreq = cos (freq);
  503. sfreq = sin (freq);
  504. cosfreq8 = cos (freq * 0.125);
  505. sinfreq8 = sin (freq * 0.125);
  506. c = cosfreq8;
  507. s = sinfreq8;
  509. for (i = 0; i < (N >> 2); i++) {
  510. /* calculate real and imaginary parts of g(n) or G(p) */
  511. tempr = -data[2 * i];
  512. tempi = data[(N >> 1) - 1 - 2 * i];
  514. /* calculate pre-twiddled FFT input */
  515. xr[i] = tempr * c - tempi * s;
  516. xi[i] = tempi * c + tempr * s;
  518. /* use recurrence to prepare cosine and sine for next value of i */
  519. cold = c;
  520. c = c * cfreq - s * sfreq;
  521. s = s * cfreq + cold * sfreq;
  522. }
  524.     /* Perform in-place complex IFFT of length N/4 */
  525. switch (N) {
  526.     case 256:
  527. srifft(xr, xi, 6);
  528. break;
  529.     case 2048:
  530. srifft(xr, xi, 9);
  531. }
  533.     /* prepare for recurrence relations in post-twiddle */
  534.     c = cosfreq8;
  535.     s = sinfreq8;
  537.     /* post-twiddle FFT output and then get output data */
  538.     for (i = 0; i < (N >> 2); i++) {
  540. /* get post-twiddled FFT output  */
  541. tempr = fac * (xr[i] * c - xi[i] * s);
  542. tempi = fac * (xi[i] * c + xr[i] * s);
  544. /* fill in output values */
  545. data [(N >> 1) + (N >> 2) - 1 - 2 * i] = tempr;
  546. if (i < (N >> 3))
  547. data [(N >> 1) + (N >> 2) + 2 * i] = tempr;
  548. else
  549. data [2 * i - (N >> 2)] = -tempr;
  551. data [(N >> 2) + 2 * i] = tempi;
  552. if (i < (N >> 3))
  553. data [(N >> 2) - 1 - 2 * i] = -tempi;
  554. else
  555. data [(N >> 2) + N - 1 - 2*i] = tempi;
  557. /* use recurrence to prepare cosine and sine for next value of i */
  558. cold = c;
  559. c = c * cfreq - s * sfreq;
  560. s = s * cfreq + cold * sfreq;
  561. }
  563. if (xr) FreeMemory(xr);
  564. if (xi) FreeMemory(xi);
  565. }