最近在做传感信号处理遇到了循环相关问题，快速循环相关算法，先将两组信号FFT变换到频域，进行复数乘积，再IFFT逆变回来即可。傅立叶变换是信号处理的基本模块，这里将这段代码拿出来分享一下，希望对同行有所用。
Dit2FFT.c文件
/********************************************************************************
*File name : dit2fft.c
*Deion : FFT calculation. Decimation in time FFT, radix-2.
*
*Current version: v2.3 -- get a simple version, remove unnecessary code
*History version:
* v2.2 -- modify the sin/cos table data, and the static
* v2.1 -- modify the name
* v2.0 -- separate the bit-reverse block
* v1.9 -- rewrite the bit-reverse block
* v1.8 -- rewrite the twiddle factor generation, and test fft,ifft
* v1.7 -- fix the twiddle factor bug, and test the four
* combinations
* v1.6 -- adjust the DIT2_FLOAT to DIT2_FLOAT32
* v1.5 -- change DIT2_UINT16 to DIT2_UINT, and extend its scope
* v1.4 -- change the name of s
* v1.3 -- add the sin/cos table
* v1.2 -- add ifft
* v1.1 -- add bit-reverse in-place
* v1.0 -- achieve the basic of FFT
*Author : Wiracle
*Date :2008-10-23
********************************************************************************/
/*inlcude the head file*/
/* the NULL define needed*/
//#include <stdlib.h>
/*
* caculate the magnitude and phase
* -- sqrt()
* -- atan()
*/
#include <math.h>
#include "Alg_Dit2FFT.h"
/*the sin(x)/cos(x) table*/
DIT2_FLOAT sin_table[21] = {
0.0, // sin(-2*pi/(2^0))
0.0, // sin(-2*pi/(2^1))
-1.0, // sin(-2*pi/(2^2))
-0.70710678118655, // sin(-2*pi/(2^3))
-0.38268343236509, // sin(-2*pi/(2^4))
-0.19509032201613, // sin(-2*pi/(2^5))
-0.09801714032956, // sin(-2*pi/(2^6))
-0.04906767432742, // sin(-2*pi/(2^7))
-0.02454122852291, // sin(-2*pi/(2^8))
-0.01227153828572, // sin(-2*pi/(2^9))
-0.00613588464915, // sin(-2*pi/(2^10))
-0.00306795676297, // sin(-2*pi/(2^11))
-0.00153398018628, // sin(-2*pi/(2^12))
-0.00076699031874, // sin(-2*pi/(2^13))
-0.00038349518757, // sin(-2*pi/(2^14))
-0.00019174759731, // sin(-2*pi/(2^15))
-0.00009587379910, // sin(-2*pi/(2^16))
-0.00004793689960, // sin(-2*pi/(2^17))
-0.00002396844981, // sin(-2*pi/(2^18))
-0.00001198422491, // sin(-2*pi/(2^19))
-0.00000599211245 // sin(-2*pi/(2^20))
};
DIT2_FLOAT cos_table[21] = {
1.0, //cos(-2*pi/(2^0))
-1.0, //cos(-2*pi/(2^1))
0.0, //cos(-2*pi/(2^2))
0.70710678118655, //cos(-2*pi/(2^3))
0.92387953251129, //cos(-2*pi/(2^4))
0.98078528040323, //cos(-2*pi/(2^5))
0.99518472667220, //cos(-2*pi/(2^6))
0.99879545620517, //cos(-2*pi/(2^7))
0.99969881869620, //cos(-2*pi/(2^8))
0.99992470183914, //cos(-2*pi/(2^9))
0.99998117528260, //cos(-2*pi/(2^10))
0.99999529380958, //cos(-2*pi/(2^11))
0.99999882345170, //cos(-2*pi/(2^12))
0.99999970586288, //cos(-2*pi/(2^13))
0.99999992646572, //cos(-2*pi/(2^14))
0.99999998161643, //cos(-2*pi/(2^15))
0.99999999540411, //cos(-2*pi/(2^16))
0.99999999885103, //cos(-2*pi/(2^17))
0.99999999971276, //cos(-2*pi/(2^18))
0.99999999992819, //cos(-2*pi/(2^19))
0.99999999998205 //cos(-2*pi/(2^20))
};
static DIT2_UINT DIT2_Idx(DIT2_UINT in);
/********************************************************************************
* name : DIT2_Rvs
*Input : DIT2_FFT pointer
*Deion : Bitreverse the input data, including real and imag part
* typedef struct DIT2_FFT_STRUCT
* {
* DIT2_FLOAT *real_in; --input and output buffer
* DIT2_FLOAT *imag_in; --input and output buffer(can be NULL)
* DIT2_FLOAT *real_out; --unused
* DIT2_FLOAT *imag_out; --unused
* DIT2_FLOAT *mag; --unused
* DIT2_FLOAT *phase; --unused
* DIT2_UINT size; --the FFT size, must be 2^N
* DIT2_UINT stage; --unused
* DIT2_UINT transform; --unused
* }DIT2_FFT;
********************************************************************************/
#pragma CODE_SECTION(DIT2_Rvs,"ramfuncs");
void DIT2_Rvs(DIT2_FFT *fft)
{
DIT2_UINT i,j,k,n;
DIT2_FLOAT tw_real,tw_imag;
/*para check*/
/*
if ((fft->real_in == NULL) || (fft->imag_in == NULL))
{
while (1);
}
*/
/*
* 1.loop along the array -- i
* 2.find the bit-reversed index -- k
* 3.check if i < k, if so, swap the X and X[k]
*/
for (i = 0; i < fft->size; i++)
{
k = 0;
n = i;
for (j = 0; j < fft->stage; j++)
{
k = (k<<1) | (n&1);
n >>= 1;
}
if (i < k)
{
tw_real = fft->real_in;
tw_imag = fft->imag_in;
fft->real_in = fft->real_in[k];
fft->imag_in = fft->imag_in[k];
fft->real_in[k] = tw_real;
fft->imag_in[k] = tw_imag;
}
}
}
/********************************************************************************
* name : DIT2_Fft
*Input : DIT2_FFT *fft
*Deion : used when uncommenting DIT2_BIT_INVERSE_IN_PLACE which achives
* the bit-reversed in-place. As the following, the real part of
* the inputsequence buffer is pointed by real_in, where the
* output--magnitude info is also in. The imaginary part of the
* input sequence buffer is pointed by image_in, where the output--
* phase info is alse in.
* typedef struct DIT2_FFT_STRUCT
* {
* DIT2_FLOAT *real_in; --input and output buffer
* DIT2_FLOAT *imag_in; --input and output buffer
* DIT2_FLOAT *real_out; --unused
* DIT2_FLOAT *imag_out; --unused
* DIT2_FLOAT *mag; --unused
* DIT2_FLOAT *phase; --unused
* DIT2_UINT size; --the FFT size, must be 2^N
* DIT2_UINT stage; --N
* DIT2_UINT transform; --DIT2_FORWARD_TRANSFORM or DIT2_INVERSE_TRANSFORM
* }DIT2_FFT;
********************************************************************************/
#pragma CODE_SECTION(DIT2_Fft,"ramfuncs");
void DIT2_Fft(DIT2_FFT *fft)
{
DIT2_UINT i,j,k,n;
DIT2_UINT block_size,block_end;
DIT2_FLOAT cos1,sin1,cos2,sin2;
DIT2_FLOAT tw_real,tw_imag; //temperary t-w
DIT2_INT16 ang_factor;
/*para check*/
/*
if ((fft->imag_in == NULL) || (fft->real_in == NULL))
{
while (1);
}
*/
/*transform direction*/
switch (fft->transform)
{
case DIT2_FORWARD_TRANSFORM:
ang_factor = DIT2_FORWARD_TRANSFORM;
break;
case DIT2_INVERSE_TRANSFORM:
ang_factor = DIT2_INVERSE_TRANSFORM;
break;
default:
break;
}
/*fft*/
block_end = 1;
/*
* The external loop:
* find the cos(t-w)/sin(t-w) -- t-w = -2*pi/N
*/
for (block_size = 2; block_size <= fft->size; block_size <<= 1)
{
i = DIT2_Idx(block_size);
cos1 = cos_table;
sin1 = sin_table;
if (ang_factor == DIT2_INVERSE_TRANSFORM)
{
sin1 = -sin1;
}
/*
* prepare the first t-w
*/
for (i = 0; i < fft->size; i += block_size)
{
cos2 = 1;//cos(0)
sin2 = 0;//sin(0)
for (j = i, n = 0; n < block_end; j++, n++)
{
k = j + block_end;
tw_real = fft->real_in[k]*cos2 - fft->imag_in[k]*sin2;
tw_imag = fft->real_in[k]*sin2 + fft->imag_in[k]*cos2;
/*lower elem*/
fft->real_in[k] = fft->real_in[j] - tw_real;
fft->imag_in[k] = fft->imag_in[j] - tw_imag;
/*uper elem*/
fft->real_in[j] = fft->real_in[j] + tw_real;
fft->imag_in[j] = fft->imag_in[j] + tw_imag;
/*adjust cos(t-w)/sin(t-w)*/
tw_real = cos2*cos1 - sin2*sin1;//cos(n+1) = cos(n)cos(1) - sin(n)sin(1)
tw_imag = cos2*sin1 + sin2*cos1;//sin(n+1) = sin(n)cos(1) + sin(1)cos(n)
cos2 = tw_real;//cos2 = cos(n+1)
sin2 = tw_imag;//sin2 = sin(n+1)
}
}
block_end = block_size;
}
/*inverse transform*/
if (fft->transform == DIT2_INVERSE_TRANSFORM)
{
for (i = 0; i < fft->size; i++)
{
fft->real_in /= fft->size;
fft->imag_in /= fft->size;
}
}
return;
}
/********************************************************************************
* name : DIT2_Mag
*Input : DIT2_FFT pointer
*Deion : calculate the magnitude and phase info of the output sequence,
* used when uncommenting the MACROA--
* #define DIT2_BIT_INVERSE_IN_PLACE
*
* typedef struct DIT2_FFT_STRUCT
* {
* DIT2_FLOAT *real_in; --output magnitude info
* DIT2_FLOAT *imag_in; --output phase info
* DIT2_FLOAT *real_out; --unused
* DIT2_FLOAT *imag_out; --unused
* DIT2_FLOAT *mag; --unused
* DIT2_FLOAT *phase; --unused
* DIT2_UINT size; --the FFT size, must be 2^N
* DIT2_UINT stage; --unused
* DIT2_UINT transform; --unused
* }DIT2_FFT;
********************************************************************************/
#pragma CODE_SECTION(DIT2_Mag,"ramfuncs");
void DIT2_Mag(DIT2_FFT *fft)
{
DIT2_UINT i;
DIT2_FLOAT mag,rad;
/*para check*/
/*
if ((fft->real_in == NULL) ||(fft->imag_in == NULL))
{
while (1);
}*/
for (i = 0; i < fft->size; i++)
{
mag = sqrt(fft->real_in * fft->real_in + fft->imag_in * fft->imag_in);
rad = atan2(fft->imag_in,fft->real_in);
fft->real_in = mag*2.0/fft->size;
fft->imag_in = rad;
}
fft->real_in[0] /= 2.0;
return;
}
/********************************************************************************
* name : DIT2_Idx
*Input : DIT2_UINT in -- the blocksize value
*Output : the index in sin_table/cos_table
*Deion : calculate the index in sin_table/cos_table to lookup the table
* according the current blocksize. Only used when uncomment the
* MACRO--#define DIT2_SIN_COS_TABALE
********************************************************************************/
#pragma CODE_SECTION(DIT2_Idx,"ramfuncs");
static DIT2_UINT DIT2_Idx(DIT2_UINT in)
{
DIT2_UINT i;
for (i=0; ; i++)
{
if (in & (1 << i))
{
return(i);
}
}
}
Dit2FFT.h文件
/********************************************************************************
*File name : Alg_Dit2FFT.h
*Deion : FFT calculation. Decimation in time FFT, radix-2.
*
*Current version: v2.3 -- get a simple version, remove unnecessary code
*History version:
* v2.2 -- modify the sin/cos table data, and the static
* v2.1 -- modify the name
* v2.0 -- separate the bit-reverse block
* v1.9 -- rewrite the bit-reverse block
* v1.8 -- rewrite the twiddle factor generation, and test fft,ifft
* v1.7 -- fix the twiddle factor bug, and test the four
* combinations
* v1.6 -- adjust the DIT2_FLOAT to DIT2_FLOAT32
* v1.5 -- change DIT2_UINT16 to DIT2_UINT, and extend its scope
* v1.4 -- change the name of s
* v1.3 -- add the sin/cos table
* v1.2 -- add ifft
* v1.1 -- add bit-reverse in-place
* v1.0 -- achieve the basic of FFT
*Author : Wiracle
*Date :2008-10-23
********************************************************************************/
#ifndef _ALG_DIT2FFT_H_
#define _ALG_DIT2FFT_H_
//#include "App_TypeDefDataStruct.h"
typedef unsigned int DIT2_UINT16;
typedef int DIT2_INT16;
typedef unsigned long DIT2_UINT32;
typedef long DIT2_INT32;
typedef long long DIT2_INT64;
typedef unsigned long long DIT2_UINT64;
typedef float DIT2_FLOAT32;
typedef long double DIT2_FLOAT64;
/*config the type used */
typedef DIT2_FLOAT32 DIT2_FLOAT;
typedef DIT2_UINT32 DIT2_UINT;
typedef DIT2_INT32 DIT2_INT;
//FFT运算数据结构
typedef struct st_DIT2_FFT
{
DIT2_FLOAT32 *real_in; //input and output buffer
DIT2_FLOAT32 *imag_in; //input and output buffer
DIT2_UINT16 size; //FFT size ,must be 2^N
DIT2_UINT16 stage; //蝶形次数
DIT2_UINT16 transform; //0: FFT运算 1：IFFT运算
}DIT2_FFT;
/*choose the FFT direction*/
#define DIT2_FORWARD_TRANSFORM 0U
#define DIT2_INVERSE_TRANSFORM 1U
//#include "..\Sources_Application\App_TypeDefDataStruct.h"
extern DIT2_FLOAT sin_table[];
extern DIT2_FLOAT cos_table[];
/********************************************************************************
* name : DIT2_Rvs
*Input : DIT2_FFT pointer
*Deion : Bitreverse the input data, including real and imag part
********************************************************************************/
extern void DIT2_Rvs(DIT2_FFT *fft);
/********************************************************************************
* name : DIT2_Fft
*Input : DIT2_FFT *fft
*Deion : Calculate the FFT or IFFT
********************************************************************************/
extern void DIT2_Fft(DIT2_FFT *fft);
/********************************************************************************
* name : DIT2_Mag
*Input : DIT2_FFT pointer
*Deion : calculate the magnitude and phase info of the output sequence,
********************************************************************************/
extern void DIT2_Mag(DIT2_FFT *fft);
#endif // end of #ifndef _ALG_DIT2FFT_H_
app_main.c 文件
// TI File $Revision: /main/14 $
// Checkin $Date: April 21, 2008 15:41:07 $
//###########################################################################
#include "DSP28x_Project.h" // Device Headerfile and Examples Include File
#include "SoftTimer.h"
#include "Alg_Dit2FFT.h"
#include "stdlib.h"
#define PI 3.1415926
// Prototype statements for s found within this file.
interrupt void cpu_timer1_isr(void);
SOFT_TIMER_LINK timer_1;
SOFT_TIMER_LINK timer_2;
#pragma DATA_SECTION(sig_real,"ZONE7DATA");
float sig_real[8192];
#pragma DATA_SECTION(sig_imag,"ZONE7DATA");
float sig_imag[8192];
DIT2_FFT dit2_fft;
void init_fft(void);
void init_ifft(void);
void main(void)
{
init_fft();
DIT2_Rvs(&dit2_fft);
DIT2_Fft(&dit2_fft);
//DIT2_Mag(&dit2_fft);
init_ifft();
DIT2_Rvs(&dit2_fft);
DIT2_Fft(&dit2_fft);
for(;;)
{
}
}
void init_fft(void)
{
int i;
float omg_o;
dit2_fft.real_in = sig_real;
dit2_fft.imag_in = sig_imag;
dit2_fft.transform = 1;
dit2_fft.size = 8192;
dit2_fft.stage = 13;
omg_o = 2*PI/dit2_fft.size;
for(i=0;i<dit2_fft.size;i++)
{
sig_real = sin(2*omg_o*i) + 2*sin(50*omg_o*i);// + 4*sin(100*omg_o*i)+ 8*sin(1000*omg_o*i) ;
sig_imag = 0;
}
}
void init_ifft(void)
{
dit2_fft.real_in = sig_real;
dit2_fft.imag_in = sig_imag;
dit2_fft.transform = 0;
dit2_fft.size = 8192;
dit2_fft.stage = 13;
}
//===========================================================================
// No more.
//===========================================================================
STM32
MCU
通讯及无线技术
物联网技术
电子DIY
板卡试用
基础知识
软件与操作系统
我爱生活
小e食堂