📄 fft.h
字号:
/* fix_fft.c - Fixed-point in-place Fast Fourier Transform */
/*
All data are fixed-point short integers, in which -32768
to +32768 represent -1.0 to +1.0 respectively. Integer
arithmetic is used for speed, instead of the more natural
floating-point.
For the forward FFT (time -> freq), fixed scaling is
performed to prevent arithmetic overflow, and to map a 0dB
sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq
coefficients. The return value is always 0.
Henceforth "short" implies 16-bit word. If this is not
the case in your architecture, please replace "short"
with a type definition which *is* a 16-bit word.
所有数据均定点短整数，用-32768~32768代表-1.0至+1.0。整
形数运算加速运算，代替浮点运算。
前向FFT中，采用固定缩放以防止算术溢出，并映射零分贝正弦/
余弦波（即振幅= 32767 ）成两个6分贝频率系数。
*/
#define N_WAVE 1024 /* full length of Sinewave[] */
#define LOG2_N_WAVE 10 /* log2(N_WAVE) */
#include "reverse_table.h"
#include "SIN_WAVE.h"
inline short FIX_MPY(short a, short b)
{
int c = ((int)a * (int)b ) >> 14;
b = c & 0x01;
a = (c >> 1) + b;
return a;
}
//=========fft核心算法=================
void fix_fft1(short fr[], short fi[], short m)
{
int i, j, k, power_2, N,a,b,c;
short qr,qi,tr,ti,wr, wi;
short *r;
N = 1 << m;
////////////////洗序////////////////
switch(N)
{
case 4: { r = r_Table4; a = 1;break;}
case 8: { r = r_Table8; a = 3;break;}
case 16: { r = r_Table16; a = 11;break;}
case 32: { r = r_Table32; a = 23;break;}
case 64: { r = r_Table64; a = 55;break;}
case 128: { r = r_Table128; a = 111;break;}
case 256: { r = r_Table256; a = 239;break;}
case 512: { r = r_Table512; a = 479;break;}
case 1024: { r = r_Table1024; a = 991;break;}
}
for(i=0;i<a;i=i+2)
{
tr = fr[r[i]];
fr[r[i]] = fr[r[i+1]];
fr[r[i+1]] = tr;
ti = fi[r[i]];
fi[r[i]] = fi[r[i+1]];
fi[r[i+1]] = ti;
}
//////////////fft核心算法/////////////////
for(i=0;i<m;i++)
{
power_2 = 1 << i;
for(j=0;j<N;j=j+power_2*2)
{
for(k=0;k<power_2;k++)
{
a = k*N_WAVE/(2*power_2);
b = j + k;
c = b + power_2;
wr = Sinewave[a+N_WAVE/4];
wi = - Sinewave[a];
tr = FIX_MPY(wr,fr[c]) - FIX_MPY(wi,fi[c]);
ti = FIX_MPY(wr,fi[c]) + FIX_MPY(wi,fr[c]);
qr = fr[b];
qi = fi[b];
//short型溢出报错,-32767~32767
if((qr-tr>32767)|(qi-ti>32767)|(qr+tr>32767)|(qi+ti>32767))
{
printf("overflow !!!!\n");
break;
}
fr[c] = qr - tr;
fi[c] = qi - ti;
fr[b] = qr + tr;
fi[b] = qi + ti;
}
}
}
}
//====================================================
// fix_fftr1----未整合函数，采用sin表,无放缩
//====================================================
short fix_fftr1(short f[], short m)
{
int i, N = 1<<m, scale = 0; //N = 128
short *fr = f;
short *fi = new short[N];
for(i=0;i<N;i++)
fi[i] = 0;
fix_fft1(fr, fi, m);
for(i=0; i<N; i++)
{
f[i] = sqrt(pow(double(fr[i]),2)+pow(double(fi[i]),2));
}
return scale;
}
//====================================================
// fix_fftr2----整合函数，采用sin表,有放缩
//====================================================
void fix_fftr2(short f[], short m)
{
unsigned short i,j,k,a,b,c,power_2, N = 1<<m;
short qr,qi,tr,ti,wr, wi;
short *fr = f;
short *fi = new short[N];
for(i=0;i<N;i++)
fi[i] = 0;
short *r;
////////////////洗序////////////////
switch(N)
{
case 4: { r = r_Table4; a = 1;break;}
case 8: { r = r_Table8; a = 3;break;}
case 16: { r = r_Table16; a = 11;break;}
case 32: { r = r_Table32; a = 23;break;}
case 64: { r = r_Table64; a = 55;break;}
case 128: { r = r_Table128; a = 111;break;}
case 256: { r = r_Table256; a = 239;break;}
case 512: { r = r_Table512; a = 479;break;}
case 1024: { r = r_Table1024; a = 991;break;}
}
for(i=0;i<a;i=i+2)
{
tr = fr[r[i]];
fr[r[i]] = fr[r[i+1]];
fr[r[i+1]] = tr;
ti = fi[r[i]];
fi[r[i]] = fi[r[i+1]];
fi[r[i+1]] = ti;
}
//////////////fft核心算法/////////////////
for(i=0;i<m;i++)
{
power_2 = 1 << i;
for(j=0;j<N;j=j+power_2*2)
{
for(k=0;k<power_2;k++)
{
a = k*N_WAVE/(2*power_2);
b = j + k;
c = b + power_2;
wr = (Sinewave[a+N_WAVE/4]>>2);
wi = - (Sinewave[a]>>2);
tr = FIX_MPY(wr,fr[c]) - FIX_MPY(wi,fi[c]);
ti = FIX_MPY(wr,fi[c]) + FIX_MPY(wi,fr[c]);
qr = (fr[b]>>2);
qi = (fi[b]>>2);
//short型溢出报错,-32767~32767
if((qr-tr>32767)|(qi-ti>32767)|(qr+tr>32767)|(qi+ti>32767))
{
printf("overflow !!!!\n");
break;
}
fr[c] = qr - tr;
fi[c] = qi - ti;
fr[b] = qr + tr;
fi[b] = qi + ti;
}
}
}
///////////////////////////////
for(i=0; i<N; i++)
{
f[i] = sqrt(pow(double(fr[i]),2)+pow(double(fi[i]),2));
f[i] = 3096/(f[i] + 1);
if(f[i]>=1548)
f[i] = 0;
}
}
//====================================================
// fix_fftr3----整合函数，采用sin表,无放缩
//====================================================
void fix_fftr3(short f[], short m)
{
unsigned short i,j,k,a,b,c,power_2, N = 1<<m;
short qr,qi,tr,ti,wr, wi;
short *fr = f;
short *fi = new short[N];
for(i=0;i<N;i++)
fi[i] = 0;
short *r;
////////////////洗序////////////////
switch(N)
{
case 4: { r = r_Table4; a = 1;break;}
case 8: { r = r_Table8; a = 3;break;}
case 16: { r = r_Table16; a = 11;break;}
case 32: { r = r_Table32; a = 23;break;}
case 64: { r = r_Table64; a = 55;break;}
case 128: { r = r_Table128; a = 111;break;}
case 256: { r = r_Table256; a = 239;break;}
case 512: { r = r_Table512; a = 479;break;}
case 1024: { r = r_Table1024; a = 991;break;}
}
for(i=0;i<a;i=i+2)
{
tr = fr[r[i]];
fr[r[i]] = fr[r[i+1]];
fr[r[i+1]] = tr;
ti = fi[r[i]];
fi[r[i]] = fi[r[i+1]];
fi[r[i+1]] = ti;
}
//////////////fft核心算法/////////////////
for(i=0;i<m;i++)
{
power_2 = 1 << i;
for(j=0;j<N;j=j+power_2*2)
{
for(k=0;k<power_2;k++)
{
a = k*N_WAVE/(2*power_2);
b = j + k;
c = b + power_2;
wr = Sinewave[a+N_WAVE/4];
wi = - Sinewave[a];
tr = FIX_MPY(wr,fr[c]) - FIX_MPY(wi,fi[c]);
ti = FIX_MPY(wr,fi[c]) + FIX_MPY(wi,fr[c]);
qr = fr[b];
qi = fi[b];
//short型溢出报错,-32767~32767
if((qr-tr>32767)|(qi-ti>32767)|(qr+tr>32767)|(qi+ti>32767))
{
printf("overflow !!!!\n");
break;
}
fr[c] = qr - tr;
fi[c] = qi - ti;
fr[b] = qr + tr;
fi[b] = qi + ti;
}
}
}
///////////////////////////////
for(i=0; i<N; i++)
{
f[i] = sqrt(pow(double(fr[i]),2)+pow(double(fi[i]),2));
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -