📄 fft.c
字号:
/************************************************************
name: fft.c
description: FFT(Fast Fourier Transform)
prototype: float fft(float *x,int n)
call: fft(x,4096)
position: dsp1.mak
version: 1.0 date: 2001-12-19
author: date: 2001-12-19
check: ZhiguoTan date: 2002-07-16
approval: date:
remark: none
***********************************************************/
/*
* FFT(Fast Fourier Transform)
* It origined from TZG who programed and tested it to get a result .
* and then was modified by YSG on 2001/07/26.
* YSG packed up the file and added a test file which is fft_main.c
*/
#include "math.h"
#define N 1024
#define PI 3.1415926535897932384626433832795
#define WRCOS cos(PI/N)//0.99999971 /*0.99999970586288*/
#define WISIN sin(PI/N)//7.66990319E-4 /*7.669903187427045e-004*/
#define NWISIN -sin(PI/N)
static void bitrv2(int n, float *a)
{
int j, j1, k, k1, l, m, m2, n2;
float xr, xi;
m = n >> 2;
m2 = m << 1;
n2 = n - 2;
k = 0;
for (j = 0; j <= m2 - 4; j += 4) {
if (j < k) {
xr = a[j];
xi = a[j + 1];
a[j] = a[k];
a[j + 1] = a[k + 1];
a[k] = xr;
a[k + 1] = xi;
} else if (j > k) {
j1 = n2 - j;
k1 = n2 - k;
xr = a[j1];
xi = a[j1 + 1];
a[j1] = a[k1];
a[j1 + 1] = a[k1 + 1];
a[k1] = xr;
a[k1 + 1] = xi;
}
k1 = m2 + k;
xr = a[j + 2];
xi = a[j + 3];
a[j + 2] = a[k1];
a[j + 3] = a[k1 + 1];
a[k1] = xr;
a[k1 + 1] = xi;
l = m;
while (k >= l) {
k -= l;
l >>= 1;
}
k += l;
}
}
/*
* functions
* cdft: Complex Discrete Fourier Transform
* function prototypes
* void cdft(int, float, float, float *);
*
*
* -------- Complex DFT (Discrete Fourier Transform) --------
* [definition]
* <case1>
* X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
* <case2>
* X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
* (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
* [usage]
* <case1>
* cdft(2*n, cos(M_PI/n), sin(M_PI/n), a);
* <case2>
* cdft(2*n, cos(M_PI/n), -sin(M_PI/n), a);
* [parameters]
* 2*n :data length (int)
* n >= 1, n = power of 2
* a[0...2*n-1] :input/output data (float *)
* input data
* a[2*j] = Re(x[j]), a[2*j+1] = Im(x[j]), 0<=j<n
* output data
* a[2*k] = Re(X[k]), a[2*k+1] = Im(X[k]), 0<=k<n
* [remark]
* Inverse of
* cdft(2*n, cos(M_PI/n), -sin(M_PI/n), a);
* is
* cdft(2*n, cos(M_PI/n), sin(M_PI/n), a);
* for (j = 0; j <= 2 * n - 1; j++) {
* a[j] *= 1.0 / n;
* }
*
*/
static void cdft(int n, float wr, float wi, float *a)
{
int i, j, k, l, m;
float wkr, wki, wdr, wdi, ss, xr, xi;
float t1;
m = n;
while (m > 4) {
l = m >> 1;
wkr = 1;
wki = 0;
t1=wi*wi;t1+=t1;
wdr = 1 - t1;
wdi = wi * wr; wdi +=wdi;
ss = wdi+wdi;
wr = wdr;
wi = wdi;
for (j = 0; j <= n - m; j += m) {
i = j + l;
xr = a[j] - a[i];
xi = a[j + 1] - a[i + 1];
a[j] += a[i];
a[j + 1] += a[i + 1];
a[i] = xr;
a[i + 1] = xi;
xr = a[j + 2] - a[i + 2];
xi = a[j + 3] - a[i + 3];
a[j + 2] += a[i + 2];
a[j + 3] += a[i + 3];
a[i + 2] = wdr * xr - wdi * xi;
a[i + 3] = wdr * xi + wdi * xr;
}
for (k = 4; k <= l - 4; k += 4) {
wkr -= ss * wdi;
wki += ss * wdr;
wdr -= ss * wki;
wdi += ss * wkr;
for (j = k; j <= n - m + k; j += m) {
i = j + l;
xr = a[j] - a[i];
xi = a[j + 1] - a[i + 1];
a[j] += a[i];
a[j + 1] += a[i + 1];
a[i] = wkr * xr - wki * xi;
a[i + 1] = wkr * xi + wki * xr;
xr = a[j + 2] - a[i + 2];
xi = a[j + 3] - a[i + 3];
a[j + 2] += a[i + 2];
a[j + 3] += a[i + 3];
a[i + 2] = wdr * xr - wdi * xi;
a[i + 3] = wdr * xi + wdi * xr;
}
}
m = l;
}
if (m > 2) {
for (j = 0; j <= n - 4; j += 4) {
xr = a[j] - a[j + 2];
xi = a[j + 1] - a[j + 3];
a[j] += a[j + 2];
a[j + 1] += a[j + 3];
a[j + 2] = xr;
a[j + 3] = xi;
}
}
if (n > 4) {
bitrv2(n, a);
}
}
/*
* functions
* rdft: Real Discrete Fourier Transform
* function prototypes
* void rdft(int, float, float, float *);
*
*
* -------- Real DFT / Inverse of Real DFT --------
* [definition]
* <case1> RDFT
* R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
* I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
* <case2> IRDFT (excluding scale)
* a[k] = R[0]/2 + R[n/2]/2 +
* sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
* sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
* [usage]
* <case1>
* rdft(n, cos(M_PI/n), sin(M_PI/n), a);
* <case2>
* rdft(n, cos(M_PI/n), -sin(M_PI/n), a);
* [parameters]
* n :data length (int)
* n >= 2, n = power of 2
* a[0...n-1] :input/output data (float *)
* <case1>
* output data
* a[2*k] = R[k], 0<=k<n/2
* a[2*k+1] = I[k], 0<k<n/2
* a[1] = R[n/2]
* <case2>
* input data
* a[2*j] = R[j], 0<=j<n/2
* a[2*j+1] = I[j], 0<j<n/2
* a[1] = R[n/2]
* [remark]
* Inverse of
* rdft(n, cos(M_PI/n), sin(M_PI/n), a);
* is
* rdft(n, cos(M_PI/n), -sin(M_PI/n), a);
* for (j = 0; j <= n - 1; j++) {
* a[j] *= 2.0 / n;
* }
*
*/
static void rdft(int n, float wr, float wi, float *a)
{
int j, k;
float wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
if (n > 4) {
wkr = 0;
wki = 0;
wdr = wi * wi;
wdi = wi * wr;
ss = wdi+wdi;ss+=ss;
wr = 1 - wdr-wdr;
wi = wdi+wdi;
if (wi >= 0) {
cdft(n, wr, wi, a);
xi = a[0] - a[1];
a[0] += a[1];
a[1] = xi;
}
for (k = (n >> 1) - 4; k >= 4; k -= 4) {
j = n - k;
xr = a[k + 2] - a[j - 2];
xi = a[k + 3] + a[j - 1];
yr = wdr * xr - wdi * xi;
yi = wdr * xi + wdi * xr;
a[k + 2] -= yr;
a[k + 3] -= yi;
a[j - 2] += yr;
a[j - 1] -= yi;
wkr += ss * wdi;
wki += ss * (0.5 - wdr);
xr = a[k] - a[j];
xi = a[k + 1] + a[j + 1];
yr = wkr * xr - wki * xi;
yi = wkr * xi + wki * xr;
a[k] -= yr;
a[k + 1] -= yi;
a[j] += yr;
a[j + 1] -= yi;
wdr += ss * wki;
wdi += ss * (0.5 - wkr);
}
j = n - 2;
xr = a[2] - a[j];
xi = a[3] + a[j + 1];
yr = wdr * xr - wdi * xi;
yi = wdr * xi + wdi * xr;
a[2] -= yr;
a[3] -= yi;
a[j] += yr;
a[j + 1] -= yi;
if (wi < 0) {
a[1] = 0.5 * (a[0] - a[1]);
a[0] -= a[1];
cdft(n, wr, wi, a);
}
} else {
if (wi < 0) {
a[1] = 0.5 * (a[0] - a[1]);
a[0] -= a[1];
}
if (n > 2) {
xr = a[0] - a[2];
xi = a[1] - a[3];
a[0] += a[2];
a[1] += a[3];
a[2] = xr;
a[3] = xi;
}
if (wi >= 0) {
xi = a[0] - a[1];
a[0] += a[1];
a[1] = xi;
}
}
/* compute the sqrt(x^2+y^2) *//*取绝对值*/
//a[0] = (a[0] > 0 ? a[0]:(-a[0]));
// for(j = 0; j < n/2; j++){
// temp = a[2*j] * a[2*j] + a[2*j+1] * a[2*j+1];
// a[j] = sqrt(temp);
//}
}
/************************************************************
name: fft.c
*************************************************************
description: Hamming window and remove DC and reduce amplitude
*************************************************************
prototype: void hamming_win_proc0(float *x,int n,float tx)
*************************************************************
reference: none
*************************************************************
revision history:
date content by version
03MAR02 create tom yue 3.01
25JUN02 comment tom yue 3.02
*************************************************************
algorithm: tx*(0.54-0.46*cos(2*pi/n))*x
************************************************************/
/*
void hamming_win_proc0(
float *x, // input/output the input data for window process
int n, // input the point number of fft
float tx)
{
int i;
float omg=2*3.1415926/n,t0,t1,t3,omgx;
t0=tx*0.54;
t1=tx*0.46;
t3=0;
for(i=0;i<4096;i++) t3+=x[i];
t3/=4096.0;
for(i=0;i<4096;i++) x[i]-=t3;
omgx=0;//tx=t3;
for (i = 0;i < n;i++){
t3=cos(omg*i);
omgx+=omg;
t3 *=t1;
t3=t0-t3;
x[i]*=t3;
}
t3=0;
for(i=0;i<4096;i++) t3+=x[i];
t3/=4096.0;
for(i=0;i<4096;i++) x[i]-=t3;
}
*/
/*
* Name : fft
*/
void fft(
float *x, /* input/output the input data for computing fft */
int n, /* input the point number of fft */
int inverse
)
{ int j; double temp;
int t=n/2;
/* Then compute the FFT */
if(!inverse){
rdft(n,WRCOS,WISIN,x);
/* compute the sqrt(x^2+y^2) 取模*/
x[0] = (x[0] > 0 ? x[0]:(-x[0]));
for(j = 0; j < n/2; j++){
temp = x[2*j] * x[2*j] + x[2*j+1] * x[2*j+1];
x[j] = sqrt(temp);
} /*取绝对值*/
}else{
rdft(n,WRCOS,NWISIN,x);
for(j=0;j<n;j++)
x[j]/=t;
}
}
////////////////////////file end////////////////////////
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -