【创龙TMS320C665x开发板试用】FFT快速傅里叶变换与IFFT逆变换

只看该作者 · 2016-10-14 17:23

利用TI DSPLIB很容易实现FFT 与IFFT

/****************************************************************************/
/* */
/* 快速傅里叶变换 / 快速傅里叶逆变换（实数） */
/* */
/* 2015年07月21日 */
/* */
/****************************************************************************/
#include <stdio.h> // C 语言标准输入输出函数库
#include <math.h> // C 数学函数库
#include "mathlib.h" // DSP 数学函数库
#include "dsplib.h" // DSP 函数库
/****************************************************************************/
/* */
/* 宏定义 */
/* */
/****************************************************************************/
// 软件断点
#define SW_BREAKPOINT asm(" SWBP 0 ");
// 快速傅里叶变换
// π 及浮点数极小值
#define PI 3.141592654
#define F_TOL (1e-06)
/****************************************************************************/
/* */
/* 全局变量 */
/* */
/****************************************************************************/
// 快速傅里叶变换测试
// 测试快速傅里叶变换点数
// 注意:TI DSP库最大支持一次性计算 256K 个实数点的 FFT
#define N 1024
// 采样频率
#define Fs 1000.0
// 信号
#pragma DATA_ALIGN(Input, 8);
float Input[2 * N];
// 信号副本
float InputOrig[2 * N];
// FFT 输出
#pragma DATA_ALIGN(FFT_Out, 8);
float FFT_Out[2 * N];
// IFFT 输出
#pragma DATA_ALIGN(IFFT_Out, 8);
float IFFT_Out[2 * N];
// 旋转因子
#pragma DATA_ALIGN(W, 8);
float W[2 * N];
// 二进制位翻转
#pragma DATA_ALIGN (brev, 8);
unsigned char brev[64] =
{
0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38,
0x04, 0x24, 0x14, 0x34, 0x0c, 0x2c, 0x1c, 0x3c,
0x02, 0x22, 0x12, 0x32, 0x0a, 0x2a, 0x1a, 0x3a,
0x06, 0x26, 0x16, 0x36, 0x0e, 0x2e, 0x1e, 0x3e,
0x01, 0x21, 0x11, 0x31, 0x09, 0x29, 0x19, 0x39,
0x05, 0x25, 0x15, 0x35, 0x0d, 0x2d, 0x1d, 0x3d,
0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b,
0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f
};
/****************************************************************************/
/* */
/* 函数声明 */
/* */
/****************************************************************************/
/****************************************************************************/
/* */
/* 快速傅里叶变换测试 */
/* */
/****************************************************************************/
// 产生旋转因子
void tw_gen(float *w, int n)
{
int i, j, k;
for(i = 1, k = 0; i < n >> 2; i++)
{
w[k ] = sin(2 * PI * i / n);
w[k + 1] = cos(2 * PI * i / n);
k += 2;
}
for(j = 1; j <= n >> 3; j = j << 2)
{
for(i = 0; i < n >> 3; i += j)
{
w[k] = (float)sin( 4 * PI * i / n);
w[k + 1] = (float)cos( 4 * PI * i / n);
w[k + 2] = (float)sin( 8 * PI * i / n);
w[k + 3] = (float)cos( 8 * PI * i / n);
w[k + 4] = (float)sin(12 * PI * i / n);
w[k + 5] = (float)cos(12 * PI * i / n);
k += 6;
}
}
}
void tw_geni(float *w, int n)
{
int i, j, k;
for(i = 1, k = 0; i < n >> 2; i++)
{
w[k ] = sin (2 * PI * i / n);
w[k + 1] = -cos (2 * PI * i / n);
k += 2;
}
for(j = 1; j <= n >> 3; j = j << 2)
{
for(i = 0; i < n >> 3; i += j)
{
w[k] = (float) -sin ( 4 * PI * i / n);
w[k + 1] = (float) cos ( 4 * PI * i / n);
w[k + 2] = (float) -sin ( 8 * PI * i / n);
w[k + 3] = (float) cos ( 8 * PI * i / n);
w[k + 4] = (float) -sin (12 * PI * i / n);
w[k + 5] = (float) cos (12 * PI * i / n);
k += 6;
}
}
}
/****************************************************************************/
/* */
/* 主函数 */
/* */
/****************************************************************************/
int main(void)
{
// 产生待测试信号
unsigned int i;
for (i = 0; i < N; i++)
{
Input = 5 * sin(2 * PI * 150 * (i / Fs)) + 15 * sin(2 * PI * 350 * (i / Fs));
}
// 基
unsigned char rad;
if(N == 64 || N == 256 || N == 1024 || N == 4096 || N == 16384 || N == 65536)
{
rad = 4;
}
else if(N == 32 || N == 128 || N == 512 || N == 2048 || N == 8192 || N == 32768)
{
rad = 2;
}
else
{
printf ("不支持计算 %d 点快速傅里叶变换！\n", N);
while(1);
}
// 保留一份输入信号副本
memcpy(InputOrig, Input, 2 * N * sizeof(float));
// 产生旋转因子
tw_gen(W, N);
// FFT 计算
DSPF_sp_fftSPxSP_r2c(N, Input, W, FFT_Out, brev, rad, 0, N);
// 计算振幅
float mo[2 * N];
for(i = 0; i < N; i++)
{
if(i == 0)
{
mo = mo / N;
}
else
{
mo = sqrtsp(powsp(FFT_Out[2 * i], 2) + powsp(FFT_Out[ 2 * i + 1], 2));
mo = mo * 2 / N;
}
}
// 产生旋转因子
tw_geni(W, N);
// IFFT 计算
DSPF_sp_ifftSPxSP_c2r(N, FFT_Out, W, IFFT_Out, brev, rad, 0, N);
for(i = 0; i < N; i++)
{
if(abs(IFFT_Out - InputOrig) > F_TOL)
{
printf("失败 %4d 原始数据 %10.5f 逆变换 %10.5f\r\n", i, InputOrig, IFFT_Out);
}
}
// 断点
SW_BREAKPOINT;
}

复制代码

130110rso7743wggzt9fk9.png.thumb.jpg.png (28.97 KB )

下载附件

2016-10-14 17:20 上传

130112ykzwwo87uhgaw8kp.jpg.thumb.jpg (35.54 KB )

下载附件

2016-10-14 17:20 上传

Matlab代码
fs=1000;                                  % 采样频率
N=1024;                                     % 数据点数
n=0:N-1;t=n/fs;                            % 时间序列
x=5*sin(2*pi*150*t)+15*sin(2*pi*350*t);    % 信号
subplot(2,1,1),plot(t,x);                   % 绘制时域信号图
y=fft(x,N);                               % 对信号进行快速傅里叶变换
mag=abs(y);mag=mag*2/N;                   % 振幅
f=n*fs/N;                                  % 频率序列
subplot(2,1,2),plot(f,mag);                % 绘制振幅图
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
grid on;                                  % 绘制振幅图
复制代码

130113yz7ro616e2l8yr53.png.thumb.jpg.png (33.24 KB )

下载附件

2016-10-14 17:20 上传

本资料来源于网络，直供学习！

只看该作者 · 2016-10-14 22:22

有没有1024的FFT实现的

只看该作者 · 2016-10-14 22:23

1024fft的cmd怎么分配？

只看该作者 · 2017-3-4 14:48

请问二维fft处理图像的代码怎么写？

只看该作者 · 2017-3-4 23:13

使用DFT要快一些。

只看该作者 · 2017-3-4 23:16

DSP做运算性能较佳。

只看该作者 · 2017-3-5 18:17

这是matlab的仿真应用么

只看该作者 · 2017-3-5 18:22

请问DFT怎么使用？

只看该作者 · 2017-3-7 17:13

看到仿真的效果很好。

只看该作者 · 2017-3-7 17:14

DSP最快做到多大的FFT？

只看该作者 · 2017-3-9 11:47

代码可执行速度怎么样

只看该作者 · 2017-3-9 11:48

比较倾向于DFT算法。

只看该作者 · 2017-3-9 21:41

matlab的程序可以直接移植吗

只看该作者 · 2017-3-9 21:43

eefas 发表于 2016-10-14 22:22
有没有1024的FFT实现的

直接修改FFT的程序不行么

【创龙TMS320C665x开发板试用】FFT快速傅里叶变换与IFFT逆变换

相关下载

相关帖子