为了在自编的一个播放器上显示频谱效果，学习了一下快速傅里叶变换，做了一个定点的FFT程序。定点FFT程序主要应用到了浮点数的定点处理和复数的运算，为此自编了一个fft_math模块。
FFT具体原理参考了《数字信号处理 第2版》(Richard G. Lyons), 在第5章节中对FFT有详尽的解释和演示。
这里需要注意的是：
1. 在Openwrt_widora编译环境下左右移位均为数学移位！也正好利用这一特点来加快计算速度。
2. 由于最大整数位的限制(设置为int64_t), 计算精度/数据范围/采样点数3者之间需要进行协调控制，不然会造成计算溢出。
3. 具体定点FFT由函数mat_egiFFFT()来实现，可看其中具体说明。
4. 考虑到不同情况下的应用，fft_math模块中分别提供了浮点型和整型的波幅计算函数。比如，对于PCM S16格式的音频数据，可用整型函数mat_uintCompAmp()来得到整型的波幅值。
5. 同样, 测试程序test_fft.c提供浮点型和整型2种采样点输入样本。注意比较不同频率其波幅值和实际计算输出的波幅值。
6. 所附程序在 Openwrt_Widora编译环境下通过。
最后应用到MP3播放器上的效果还是不错的。
----------------- 参考源代码 ---------------
fft_math.h
/*------------------------ fft_math.h ---------------------------
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License version 2 as
published by the Free Software Foundation.
Midas Zhou
-------------------------------------------------------------------*/
#ifndef __FFT_MATH_H__
#define __FFT_MATH_H__
#include <inttypes.h>
#define MATH_PI 3.1415926535897932
#define MATH_DIVEXP 11 /* An odd number!! exponent of 2, as for divisor of fixed point number */
/* EGI fixed point / complex number */
typedef struct {
int64_t num; /* divident */
int div; /* divisor, 2 exponent taken 16*/
} EGI_FVAL;
typedef struct {
EGI_FVAL real; /* real part */
EGI_FVAL imag; /* imaginery part */
} EGI_FCOMPLEX;
/* TODO: adjustable exponent value MAT_FVAL(a,exp) */
/* shitf MATH_DIVEXP bitS each time, if 15, 64-15-15-15-15-1 */
/* create a fix value */
/* (a) is float */
#define MAT_FVAL(a) ( (EGI_FVAL){ (a)*(1U<<MATH_DIVEXP), MATH_DIVEXP} )
/* (a) is INT, suppose left bit_shift is also arithmatic !!! */
#define MAT_FVAL_INT(a) ( (EGI_FVAL){ a<<MATH_DIVEXP, MATH_DIVEXP} )
/* create a complex: r--real, i--imaginery */
/* (r) is float */
#define MAT_FCPVAL(r,i) ( (EGI_FCOMPLEX){ MAT_FVAL(r), MAT_FVAL(i) } )
/* (r) is INT, suppose left bit_shift is also arithmatic !!! */
#define MAT_FCPVAL_INT(r,i) ( (EGI_FCOMPLEX){ MAT_FVAL_INT(r), MAT_FVAL(i) } )
/* complex phase angle factor
@n: m*(2*PI)/N nth phase angle
@N: Total parts of a circle(2*PI).
*/
#define MAT_CPANGLE(n, N) ( MAT_FCPVAL( cos(2.0*n*MATH_PI/N), -sin(2.0*n*MATH_PI/N) ) )
/* complex and FFT functions */
void mat_FixPrint(EGI_FVAL a);
void mat_CompPrint(EGI_FCOMPLEX a);
inline float mat_floatFix(EGI_FVAL a);
inline EGI_FVAL mat_FixAdd(EGI_FVAL a, EGI_FVAL b);
inline EGI_FVAL mat_FixSub(EGI_FVAL a, EGI_FVAL b);
inline EGI_FVAL mat_FixMult(EGI_FVAL a, EGI_FVAL b);
inline EGI_FVAL mat_FixDiv(EGI_FVAL a, EGI_FVAL b);
inline EGI_FCOMPLEX mat_CompAdd(EGI_FCOMPLEX a, EGI_FCOMPLEX b);
inline EGI_FCOMPLEX mat_CompSub(EGI_FCOMPLEX a, EGI_FCOMPLEX b);
inline EGI_FCOMPLEX mat_CompMult(EGI_FCOMPLEX a, EGI_FCOMPLEX b);
inline EGI_FCOMPLEX mat_CompDiv(EGI_FCOMPLEX a, EGI_FCOMPLEX b);
float mat_floatCompAmp( EGI_FCOMPLEX a );
unsigned int mat_uint32Log2(uint32_t np);
unsigned int mat_uintCompAmp( EGI_FCOMPLEX a );
uint64_t mat_uintCompSAmp( EGI_FCOMPLEX a );
EGI_FCOMPLEX *mat_CompFFTAng(uint16_t np);
int mat_egiFFFT( uint16_t np, const EGI_FCOMPLEX *wang,
const float *x, const int *nx, EGI_FCOMPLEX *ffx);
/* other math functions */
uint64_t mat_fp16_sqrtu32(uint32_t x);
unsigned int tm_diffus(struct timeval t_start, struct timeval t_end);
#endiffft_math.c
/*---------------------- fft_math.c --------------------------------
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License version 2 as
published by the Free Software Foundation.
Midas Zhou
--------------------------------------------------------------------*/
#include <stdio.h>
#include <sys/mman.h>
#include <string.h>
#include <math.h>
#include <inttypes.h>
#include <stdlib.h>
#include "fft_math.h"
/* ------------------------ PIXED_POINT FUNCTION ------------------------
1. After calculatoin, a Fix point value must reset its exponent 'div' to the
default value as '15'
2. Limitation analysis:
1st grade precision: 0.00003
--------------------------------------------------------------------------*/
/*---------------------------
Print function
---------------------------*/
void mat_FixPrint(EGI_FVAL a)
{
printf("Float: %.8f ", mat_floatFix(a) );
printf("[num:%"PRId64", div:%d]",a.num, a.div);
}
/*-------------------------------------------------------
Convert fixed point to float value
--------------------------------------------------------*/
inline float mat_floatFix(EGI_FVAL a)
{
// return (float)1.0*a.num/(1u<<a.div);
return 1.0*(a.num>>a.div);
}
/*-------------------------------------------------------
Addup of two fixed point value: a+b
NOTE: Divisors of two EGI_FVAL must be the same!
--------------------------------------------------------*/
inline EGI_FVAL mat_FixAdd(EGI_FVAL a, EGI_FVAL b)
{
return (EGI_FVAL){ a.num+b.num, a.div };
}
/*-------------------------------------------------------
Subtraction of two fixed point value: a+b
NOTE: Divisors of two EGI_FVAL must be the same!
--------------------------------------------------------*/
inline EGI_FVAL mat_FixSub(EGI_FVAL a, EGI_FVAL b)
{
return (EGI_FVAL){ a.num-b.num, a.div };
}
/*-------------------------------------------------------
Mutliplication of two fixed point value: a+b
NOTE: Divisors of two EGI_FVAL must be the same!
--------------------------------------------------------*/
inline EGI_FVAL mat_FixMult(EGI_FVAL a, EGI_FVAL b)
{
int64_t c;
c=a.num*b.num; /* !!!! Limit: (16num+16div)bit*(16num+16div)bit, for DIVEXP=15 */
c>>=a.div; /* Right shifting is supposed to be arithmatic */
return (EGI_FVAL){c, a.div};
}
/*-------------------------------------------------------
Division of two fixed point value: a/b
NOTE: Divisors of two EGI_FVALV must be the same!
--------------------------------------------------------*/
inline EGI_FVAL mat_FixDiv(EGI_FVAL a, EGI_FVAL b)
{
/* Min value for 0 */
if(b.num==0)
b.num=1;
return (EGI_FVAL){ (a.num<<a.div)/b.num, a.div }; /* suppose left shif is arithmatic! */
}
/*-------------------- PIXED_POINT COMPLEX FUNCTION ----------------*/
/*-------------------------------
Print function
--------------------------------*/
void mat_CompPrint(EGI_FCOMPLEX a)
{
printf(" Float: %.8f %+.8fj ", mat_floatFix(a.real), mat_floatFix(a.imag));
printf(" Real[num:%"PRId64", div:%d],Imag[num:%"PRId64", div:%d] ",
a.real.num, a.real.div, a.imag.num, a.imag.div);
}
/*-------------------------------------------------------
Addup two complex: a+b
NOTE: Divisors of two EGI_FVAL must be the same!
--------------------------------------------------------*/
inline EGI_FCOMPLEX mat_CompAdd(EGI_FCOMPLEX a, EGI_FCOMPLEX b)
{
return (EGI_FCOMPLEX){
{( a.real.num+b.real.num), a.real.div },
{( a.imag.num+b.imag.num), a.imag.div }
};
}
/*-------------------------------------------------------
Subtraction of two complex: a-b
NOTE: Divisors of two EGI_FCOMPLEX must be the same!
--------------------------------------------------------*/
inline EGI_FCOMPLEX mat_CompSub(EGI_FCOMPLEX a, EGI_FCOMPLEX b)
{
return (EGI_FCOMPLEX){
{a.real.num-b.real.num, a.real.div },
{a.imag.num-b.imag.num, a.imag.div }
};
}
/*-------------------------------------------------------
Multiplication of two complex: a*b
NOTE: Divisors of two EGI_FCOMPLEX must be the same!
--------------------------------------------------------*/
inline EGI_FCOMPLEX mat_CompMult(EGI_FCOMPLEX a, EGI_FCOMPLEX b)
{
EGI_FVAL real;
EGI_FVAL imag;
real=mat_FixSub(
mat_FixMult(a.real, b.real),
mat_FixMult(a.imag, b.imag)
);
imag=mat_FixAdd(
mat_FixMult(a.real, b.imag),
mat_FixMult(a.imag, b.real)
);
return (EGI_FCOMPLEX){real, imag};
}
/*-------------------------------------------------------
Division of two complex: a/b
NOTE: Divisors of two EGI_FCOMPLEX must be the same!
ar+jai/(br+jbi)=(ar+jai)(br-jbi)/(br^2+bi^2)
--------------------------------------------------------*/
inline EGI_FCOMPLEX mat_CompDiv(EGI_FCOMPLEX a, EGI_FCOMPLEX b)
{
EGI_FCOMPLEX d,ad;
EGI_FVAL den;
/* d=br-j(bi) */
d.real=b.real;
d.imag.num=-b.imag.num;
d.imag.div=b.imag.div;
/* a*d */
ad=mat_CompMult(a,d);
/* rb*rb+ib*ib */
den=mat_FixAdd( mat_FixMult(b.real, b.real),
mat_FixMult(b.imag, b.imag) );
/* ad/den */
ad.real=mat_FixDiv(ad.real, den); /* 0 will reset to 1.0/2^15 in FixDiv() */
ad.imag=mat_FixDiv(ad.imag, den);
return ad;
}
/*---------------------------------------------------------
work out the amplitude(modulus) of a complex, in float type.
Check also: int mat_intCompAmp( EGI_FCOMPLEX a )
Limit: a.real MAX. 2^16. (for DIVEXP=15)
---------------------------------------------------------*/
float mat_floatCompAmp( EGI_FCOMPLEX a )
{
EGI_FVAL c; /* ar^2+ai^2 */
uint64_t uamp;
/* Limit here as in mat_FixMult(): (16+16)bit*(16+16)bit */
c=mat_FixAdd( mat_FixMult(a.real, a.real),
mat_FixMult(a.imag, a.imag)
);
#if 0 /* 1. use float method sqrt */
return sqrt(mat_floatFix(c));
#else /* 2. use fixed point method sqrtu32(Max.1<<29) */
uamp=mat_fp16_sqrtu32( c.num>>(47+(15-MATH_DIVEXP)-29 +1))<<((47+(15-MATH_DIVEXP)-29)/2-MATH_DIVEXP/2);
/* if float */
return 1.0*uamp/(1<<16);
#endif
}
/*--------------------------------------
Return log2 of np
@np: >0
--------------------------------------*/
unsigned int mat_uint32Log2(uint32_t np)
{
int i;
for(i=0; i<32; i++) {
if( (np<<i) & (1<<31) ){
return 31-i;
}
}
return 0;
}
/*---------------------------------------------------------
work out the amplitude(modulus) of a complex, in INT type.
check also: float mat_floatCompAmp( EGI_FCOMPLEX a )
---------------------------------------------------------*/
unsigned int mat_uintCompAmp( EGI_FCOMPLEX a )
{
EGI_FVAL c; /* ar^2+ai^2 */
uint64_t uamp;
/* Limit here as in mat_FixMult(): (16+16)bit*(16+16)bit */
c=mat_FixAdd( mat_FixMult(a.real, a.real),
mat_FixMult(a.imag, a.imag)
);
#if 0 /* 1. use float method sqrt */
return sqrt(mat_floatFix(c));
#else /* 2. use fixed point method sqrtu32(Max.1<<29) */
uamp=mat_fp16_sqrtu32( c.num>>(47+(15-MATH_DIVEXP)-29 +1))<<((47+(15-MATH_DIVEXP)-29)/2-MATH_DIVEXP/2);
return uamp>>16;
#endif
}
/*---------------------------------------------------
Return square of the amplitude(modulus) of a complex
check also: float mat_floatCompAmp( EGI_FCOMPLEX a )
----------------------------------------------------*/
uint64_t mat_uintCompSAmp( EGI_FCOMPLEX a )
{
EGI_FVAL c; /* ar^2+ai^2 */
/* Limit here as in mat_FixMult(): (16+16)bit*(16+16)bit (for DIVEXP=15) */
c=mat_FixAdd( mat_FixMult(a.real, a.real),
mat_FixMult(a.imag, a.imag)
);
return c.num>>MATH_DIVEXP;
}
/*----------------------------------------------------
Generate complex phase angle array for FFT
@np Phanse angle numbers as per FFT point numbers
np will be normalized first to powers of 2
!!! DO NOT FORGET TO FREE IT !!!
Return:
a pointer to EGI_FCOMPLEX OK
NULL fails
----------------------------------------------------*/
EGI_FCOMPLEX *mat_CompFFTAng(uint16_t np)
{
int i;
int exp=0;
int nn;
EGI_FCOMPLEX *wang;
/* get exponent number of 2 */
for(i=0; i<16; i++) {
if( (np<<i) & (1<<15) ){
exp=16-i-1;
break;
}
}
/* reset nn to be powers of 2 */
nn=1<<exp;
/* calloc wang */
wang=calloc(nn, sizeof(EGI_FCOMPLEX));
if(wang==NULL) {
printf("%s,Fail to calloc wang!\n",__func__);
return NULL;
}
/* generate phase angle */
for(i=0;i<nn;i++) {
wang[i]=MAT_CPANGLE(i, nn);
//mat_CompPrint(wang[i]);
//printf("\n");
}
return wang;
}
/*-------------------------------------------------------------------------------------
Fixed point Fast Fourier Transform
@np: Total number of data for FFT, will be ajusted to a number of 2 powers;
np=1, result is 0!
@wang complex phase angle factor, according to normalized np.
@x: Pointer to array of input float data x[].
@nx: Pointer to araay of input INT data nx[]. NOTE: if x==NULL, use nx[] then.
@ffx: Pointer to array of FFT result ouput in form of x+jy.
Note:
0. If input x[] is NULL, then use nx[].
1. Input number of points will be ajusted to a number of powers of 2.
2. Actual amplitude is FFT_A/(NN/2), where FFT_A is ffx[] result amplitude(modulus).
3. Use real and imaginery part of ffx[] to get phase angle(relative to cosine!!!):
atan(ffx[].real/ffx[].imag),
4. !!!--- Cal overflow ---!!!
Expected amplitude: MSB. 1<<N ( 1<<(N+1)-1 )
expected sampling point: 1<<M
Taken: N+M-1=16 to incorporate with mat_floatCompAmp()
OR N+M-1=16+(15-MATH_DIVEXP) --> N+M=17+(15-MATH_DIVEXP) [ MATH_DIVEXP MUST be odd!]
Example: when MATH_DIVEXP = 15 LIMIT: Amp_Max*(np/2) is 17bits;
we may take Amp_Max: (1<<(12+1))-1; FFT point number: 1<<5:
(NOTE: 12 is the highes significant bit for amplitude )
when MATH_DIVEXP = 11 LIMIT: Amp_Max*(np/2) is 22bits;
we may take Amp_Max: (1<<(11+1))-1; FFT point number: 1<<10:
Ok, with reduced precision as you decrease MATH_DIVEXP !!!!!
5. So,you have to balance between data precision, data scope and number of sample points(np)
,which also decides Min. analysis frequency.
Return:
0 OK
<0 Fail
---------------------------------------------------------------------------------------*/
int mat_egiFFFT( uint16_t np, const EGI_FCOMPLEX *wang,
const float *x, const int *nx,EGI_FCOMPLEX *ffx)
{
int i,j,s;
int kx,kp;
int exp=0; /* 2^exp=nn */
int nn; /* number of data for FFT analysis, a number of 2 powers */
EGI_FCOMPLEX *ffodd; /* for FFT STAGE 1,3,5.. , store intermediate result */
EGI_FCOMPLEX *ffeven; /* for FFT STAGE 0,2,4.. , store intermediate result */
int *ffnin; /* normal order mapped index */
int fftmp;
/* check input data */
if( np==0 || wang==NULL || ( x==NULL && nx==NULL ) || ffx==NULL) {
printf("%s: Invalid input data!\n",__func__);
return -1;
}
/* get exponent number of 2 */
for(i=0; i<16; i++) {
if( (np<<i) & (1<<15) ){
exp=16-i-1;
break;
}
}
/* reset nn to be powers of 2 */
nn=1<<exp;
/* allocate vars !!!!! This is time consuming !!!!!! */
ffodd=calloc(nn, sizeof(EGI_FCOMPLEX));
if(ffodd==NULL) {
printf("%s: Fail to alloc ffodd[] \n",__func__);
return -1;
}
ffeven=calloc(nn, sizeof(EGI_FCOMPLEX));
if(ffeven==NULL) {
printf("%s: Fail to alloc ffeven[] \n",__func__);
free(ffodd);
return -1;
}
ffnin=calloc(nn, sizeof(int));
if(ffnin==NULL) {
printf("%s: Fail to alloc ffnin[] \n",__func__);
free(ffodd); free(ffeven);
return -1;
}
/// FFT
/* 1. map normal order index to input x[] index */
for(i=0; i<nn; i++)
{
ffnin[i]=0;
for(j=0; j<exp; j++) {
/* move i(j) to i(0), then left shift (exp-j) bits and assign to ffnin[i](exp-j) */
ffnin[i] += ((i>>j)&1) << (exp-j-1);
}
//printf("ffnin[%d]=%d\n", i, ffnin[i]);
}
/* 2. store x() to ffeven[], index as mapped according to ffnin[] */
if(x) { /* use float type x[] */
for(i=0; i<nn; i++)
ffeven[i]=MAT_FCPVAL(x[ffnin[i]], 0.0);
} else { /* use INT type nx[] */
for(i=0; i<nn; i++) {
//ffeven[i]=MAT_FCPVAL_INT(nx[ffnin[i]], 0.0);
ffeven[i].real.num=nx[ffnin[i]]<<MATH_DIVEXP; /* suppose left bit_shift is also arithmatic! */
ffeven[i].real.div=MATH_DIVEXP;
ffeven[i].imag.num=0;
ffeven[i].imag.div=MATH_DIVEXP;
}
}
/* 3. stage 2^1 -> 2^2 -> 2^3 -> 2^4....->NN point DFT */
for(s=1; s<exp+1; s++) {
for(i=0; i<nn; i++) { /* i as normal order index */
/* get coupling x order index: ffeven order index -> x order index */
kx=((i+ (1<<(s-1))) & ((1<<s)-1)) + ((i>>s)<<s); /* (i+2^(s-1)) % (2^s) + (i/2^s)*2^s) */
/* get complex phase angle index */
kp= (i<<(exp-s)) & ((1<<exp)-1); // k=(i*1)%8
if( s & 1 ) { /* odd stage */
if(i < kx )
ffodd[i]=mat_CompAdd( ffeven[i], mat_CompMult(ffeven[kx], wang[kp]) );
else /* i > kx */
ffodd[i]=mat_CompAdd( ffeven[kx], mat_CompMult(ffeven[i], wang[kp]) );
}
else { /* even stage */
if(i < kx) {
ffeven[i]=mat_CompAdd( ffodd[i], mat_CompMult(ffodd[kx], wang[kp]) );
//printf(" stage 2: ffeven[%d]=ffodd[%d]+ffodd[%d]*wang[%d]\n", i,i,kx,kp);
}
else { /* i > kx */
ffeven[i]=mat_CompAdd( ffodd[kx], mat_CompMult(ffodd[i], wang[kp]) );
//printf(" stage 2: ffeven[%d]=ffodd[%d]+ffodd[%d]*wang[%d]\n", i,kx,i,kp);
}
} /* end of odd and even stage cal. */
} /* end for(i) */
} /* end for(s) */
/* 4. pass result */
if(exp&1) { /* 4.1 result in odd array */
for(i=0; i<nn; i++) {
ffx[i]=ffodd[i];
}
} else { /* 4.2 result in even array */
for(i=0; i<nn; i++) {
ffx[i]=ffeven[i];
}
}
/* free resources */
free(ffodd);
free(ffeven);
free(ffnin);
return 0;
}
/*----------------------------------------------------------------------------------
Fix point square root for uint32_t, the result is <<16 shifted, so you need to shift
>>16 back after call to get the right answer.
Limit input x <= 2^29 //4294967296(2^32) 4.2949673*10^9
Original Source refert to:
http://read.pudn.com/downloads286/sourcecode/embedded/1291109/sqrt.c__.htm
------------------------------------------------------------------------------------*/
uint64_t mat_fp16_sqrtu32(uint32_t x)
{
uint64_t x0=x;
x0<<=32;
uint64_t x1;
uint32_t s=1;
uint64_t g0,g1;
if(x<=1)return x;
x1=x0-1;
if(x1>4294967296-1) {
s+=16;
x1>>=32;
}
if(x1>65535) {
s+=8;
x1>>=16;
}
if(x1>255) {
s+=4;
x1>>=8;
}
if(x1>15) {
s+=2;
x1>>=4;
}
if(x1>3) {
s+=1;
}
g0=1<<s;
g1=(g0+(x0>>s))>>1;
while(g1<g0) {
g0=g1;
g1=(g0+x0/g0)>>1;
}
return g0; /* after call, right shift >>16 to get right answer */
}
/*------------------------------------------------------------------
return time difference in us, as an unsigned value.
t_start: start time
t_end: end time
------------------------------------------------------------------*/
unsigned int tm_diffus(struct timeval t_start, struct timeval t_end)
{
int ds=t_end.tv_sec-t_start.tv_sec;
int dus=t_end.tv_usec-t_start.tv_usec;
int td=ds*1000000+dus;
return ( td>0 ? td : -td );
}
测试程序 test_fft.c
/*------------------------------------------------------------------
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License version 2 as
published by the Free Software Foundation.
An example to test EGI Fixed Point Fast Fourier Transform function.
Midas Zhou
------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include "fft_math.h"
#include <math.h>
#include <inttypes.h>
#define FLOAT_INPUT 1 /* 1 -- input/out in float, or 0 -- in INT */
int main(void)
{
int i,k;
struct timeval tm_start,tm_end;
/* nexp+aexp MAX 21 */
int nexp=11; /* exponent of 2, for Element number */
int aexp=10; /* exponent of 2, MSB for Max. amplitude, Max amp=2^12-1 */
uint16_t np=1<<nexp;
float famp;
unsigned int namp;
unsigned int nsamp;
float *x; /* input float */
int *nx; /* input INT */
EGI_FCOMPLEX *ffx; /* result */
EGI_FCOMPLEX *wang; /* unit phase angle factor */
/* allocate vars */
x=calloc(np, sizeof(float));
if(x==NULL)
return -1;
nx=calloc(np, sizeof(int));
if(nx==NULL)
return -1;
ffx=calloc(np, sizeof(EGI_FCOMPLEX));
if(ffx==NULL)
return -1;
/* generate NN points samples */
#if FLOAT_INPUT /* --- float type x[] --- */
for(i=0; i<np; i++) {
x[i]=10.5*cos(2.0*MATH_PI*1000*i/8000) +0.5*cos(2.0*MATH_PI*2000*i/8000+3.0*MATH_PI/4.0) \
+1.0*((1<<aexp))*cos(2.0*MATH_PI*3000*i/8000+MATH_PI/4.0); /* Amplitud Max. abt.(1<<12)+(1<<10) */
printf(" x[%d]: %f\n", i, x[i]);
}
#else /* --- INT type nx[] --- */
for(i=0; i<np; i++) {
nx[i]= (int)( 333*cos(2.0*MATH_PI*1000*i/16000) )
+(int)( 444*cos(2.0*MATH_PI*2000*i/16000+3.0*MATH_PI/4.0) )
+(int)( 555*cos(2.0*MATH_PI*3000*i/16000-1.0*MATH_PI/4.0) )
+(int)( ((1<<(aexp+1))-1)*cos(2.0*MATH_PI*5000*i/16000+MATH_PI/4.0) );
printf(" nx[%d]: %d\n", i, nx[i]);
}
#endif
/* prepare phase angle */
wang=mat_CompFFTAng(np);
k=0;
do { /* ---------------------- LOOP TEST FFT ---------------------- */
k++;
gettimeofday(&tm_start, NULL);
/* call fix_FFT */
#if FLOAT_INPUT /* float x[] */
mat_egiFFFT(np, wang, x, NULL, ffx);
#else /* int nx[] */
mat_egiFFFT(np, wang, NULL, nx, ffx);
#endif
gettimeofday(&tm_end, NULL);
/* print result */
#if 1
if( (k&(64-1))==0 )
for(i=0; i<np; i++) {
#if FLOAT_INPUT /* --- get float modulus --- */
famp=mat_floatCompAmp(ffx[i])/(np/2); /* np MUST be powers of 2 */
if( famp > 0.001 || famp < -0.001 ) {
printf("ffx[%d] ",i);
mat_CompPrint(ffx[i]);
printf(" Amp=%f ", famp);
printf("\n");
}
#else /* --- get int moduls --- */
namp=mat_uintCompAmp(ffx[i])/(np/2); /* np MUST be powers of 2 */
if( namp > 0 ) {
printf("ffx[%d] ",i);
mat_CompPrint(ffx[i]);
printf(" Amp=%d ", namp);
if(nexp>1)
nsamp=mat_uintCompSAmp(ffx[i])>>(2*(nexp-1)); // /(np/2)/(np/2) );
else
nsamp=mat_uintCompSAmp(ffx[i]);
printf(" SAmp=%d\n", nsamp);
printf("\n");
}
#endif
}
#endif
/* print results every 64 rounds */
if( (k&(64-1))==0 )
printf(" --- K=%d time cost: %d us\n", k, tm_diffus(tm_start,tm_end));
} while(1); /* ------------------ END LOOP TEST ----------------- */
/* free resources */
free(wang);
free(x);
free(nx);
free(ffx);
return 0;
}
1921
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
