STC89C52RC单片机实现FFT运算，问题求助

单片机

   如题，本人用STC89C52RC单片机实现FFT运算。我想做128点FFT，从网上找了一段纯C的代码，经过修改后移植到单片机上，但是发现当点数超过16时不出结果，点数为16及以下时，结果与matlab运算结果一致，代码反复修改了很多次，还是没解决点数多时不出结果的问题，特在此提问，希望得到大家的帮助，先表示感谢。

ps: 89C52RC具有8Kflash，512字节RAM。
   1、假如我要做32点的话，占内存的主要是定义了大小为32的结构体数组，每个结构体2个float型变量，这是32*2*4=256个字节，还有一个float型数组,大小为16，这是16*4=64字节，其余变量局部的和全局的都加起来不到60字节，所以总共是不到380字节，按380字节算。单片机运行是不是要占用一部分内存？然后剩余的内存不足380字节，然后造成不出结果，所以点数为32时不出结果是不是内存不够的原因？可能的原因还有哪些？
   2、由于要在单片机上实现FFT算法，所以点数最好为2的指数，方便运算；
   3、程序在VC++6.0上跑过，16、32、64、128、256点均试过，结果均正确（就是在单片机上16点以上都跑不动）；
   3、附件是uVision工程，可以下载调试

程序流程如下：

#include
#include "fft.h"
//main函数
void main(void)
{
u8 Freq=0; //峰值对应频率
u16 Fs=100; //采样频率
//从y=200*sin(2*pi*40*t)中取得数据点
s16 code DATA[128]={0,118 , -190 , 190 , -118 , 0 , 118 , -190 ,
190 ,-118 , 0 , 118 , -190 , 190 , -118, 0 ,
118 ,-190 , 190 , -118 , 0 , 118 , -190 , 190 ,
-118,0 ,118 , -190 , 190 , -118, 0 , 118 ,
-190 , 190 , -118 , 0 , 118 , -190 , 190 , -118 ,
0 , 118 , -190 , 190 , -118, 0 , 118 , -190 ,
190 , -118 , 0 , 118 , -190 , 190 , -118 , 0 ,
118 , -190 , 190 , -118 , 0 , 118 , -190 , 190 ,
-118 , 0 , 118 , -190 , 190 , -118 , 0 , 118 ,
-190 , 190 , -118 , 0 , 118 , -190 , 190 , -118 ,
0 , 118 , -190 , 190 , -118 , 0 , 118 , -190 ,
190 , -118 , 0 , 118 , -190 , 190 , -118 , 0 ,
118 , -190 , 190 , -118 , 0 , 118 , -190 , 190 ,
-118 , 0 , 118 , -190 , 190 , -118 , 0 , 118 ,
-190 , 190 , -118 , 0 , 118 , -190 , 190 , -118,
0 , 118 , -190, 190 , -118, 0 , 118 , -190 };
u8 ii,Max_Place=0; //用Max_Place保存频谱中最大峰值在数组中的位置
float Out_abs[NN/2]; //用于保存模值，因为频谱图关于Fs/2对称，所以只求一半
complex y[NN];
for (ii = 0; ii < NN; ++ii)
{
y[ii].real = DATA[ii];
y[ii].imag = 0;
}
fft(NN, y); //FFT运算
c_abs(y, Out_abs, NN/2); //求模
Max_Place=Mag(Out_abs); //求最大模所在位置
Freq=Max_Place*Fs/NN; //得到最大模对应频率
print(Freq);
}

fft.h

#ifndef __FFT_H__
#define __FFT_H__
typedef struct complex //复数类型
{
float real; //实部
float imag; //虚部
}complex;
#define PI 3.1415926535
#define NN 128 //采样点数,此处变更采样点数，所有用到NN的函数会自动同步更改
#define u8 unsigned char
#define s8 signed char
#define u16 unsigned int
#define s16 signed int
///////////////////////////////////////////In file fft.c
void c_plus(complex a,complex b,complex *c);//复数加
void c_mul(complex a,complex b,complex *c) ;//复数乘
void c_sub(complex a,complex b,complex *c); //复数减法
void fft(int N,complex f[]);//傅立叶变换输出也存在数组f中
void c_abs(complex f[],float out[],int n);//复数数组取模
u8 Mag(float outRes[]);
////////////////////////////////////////////
//In file display.c
void display(u8 bai,u8 shi,u8 ge);
void delayms(u8 xms);
void print(u8 num);
#endif

fft.c

#include "math.h"
#include "fft.h"
//精度0.0001弧度
void c_abs(complex f[],float out[],int n)
{
int i = 0;
float t;
for(i=0;i {
t = f.real * f.real + f.imag * f.imag;
out = sqrt(t);
}
}
void c_plus(complex a,complex b,complex *c)
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
void c_sub(complex a,complex b,complex *c)
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
void c_mul(complex a,complex b,complex *c)
{
c->real = a.real * b.real - a.imag * b.imag;
c->imag = a.real * b.imag + a.imag * b.real;
}
void Wn_i(int n,int i,complex *Wn)
{
Wn->real = cos(2*PI*i/n);
Wn->imag = -sin(2*PI*i/n);
}
//傅里叶变化
void fft(int N,complex f[])
{
complex t,wn;//中间变量
int i,j,k,m,n,l,r,M;
int la,lb,lc;
/*----计算分解的级数M=log2(N)----*/
for(i=N,M=1;(i=i/2)!=1;M++);
/*----按照倒位序重新排列原信号----*/
for(i=1,j=N/2;i<=N-2;i++)
{
if(i {
t=f[j];
f[j]=f;
f=t;
}
k=N/2;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
/*----FFT算法----*/
for(m=1;m<=M;m++)
{
la=pow(2,m); //la=2^m代表第m级每个分组所含节点数
lb=la/2; //lb代表第m级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离
/*----碟形运算----*/
for(l=1;l<=lb;l++)
{
r=(l-1)*pow(2,M-m);
for(n=l-1;n {
lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号
Wn_i(N,r,&wn);//wn=Wnr
c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算
c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
}
}
}
}
u8 Mag(float outRes[])
{
u8 i=1;
u8 max,place;
max=outRes;
place=i;
for(i=2;i {
if(outRes>max)
{
max=outRes;
place=i;
}
}
return place;
}

最后是显示函数display.c

#include
#include "fft.h"
u8 code table[]={0x3f,0x06,0x5b,0x4f,0x66,0x6d,0x7d,
0x07,0x7f,0x6f,0x77,0x7c,0x39,0x5e,0x79,0x71};
//延时函数,毫秒级
void delayms(u16 xms){
u8 j;
u16 i;
for(i=xms;i>0;i--)
for(j=110;j>0;j--);
}
//数码管显示百位十位个位数字
void display(u8 bai,u8 shi,u8 ge)
{
dula=1;
P0=table[bai];
dula=0;
P0=0xff;
wela=1;
P0=0x7e;
wela=0;
delayms(5);

dula=1;
P0=table[shi];
dula=0;
P0=0xff;
wela=1;
P0=0x7d;
wela=0;
delayms(5);
dula=1;
P0=table[ge];
dula=0;
P0=0xff;
wela=1;
P0=0x7b;
wela=0;
delayms(5);
}
//显示数字，用于调试
void print(u8 num)
{
u8 bai=0;
u8 shi=0;
u8 ge=0;
bai=num/100;
shi=num%100/10;
ge=num%10;
display(bai,shi,ge);
delayms(5000);
}

FFTWay2 uvision.zip (41.5 KB, 下载次数: 24 )