寻求快速傅立叶算法的Ｃ语言实现

air23feng · 发表于 2006-4-13 16:29:22

如题，望大虾点拨！

_yuming · 发表于 2006-4-13 16:36:34

#include "math.h"

  void kfft(pr,pi,n,k,fr,fi,l,il)

  int n,k,l,il;

  double pr[],pi[],fr[],fi[];

  { int it,m,is,i,j,nv,l0;

double p,q,s,vr,vi,poddr,poddi;

for (it=0; it<=n-1; it++)

   { m=it; is=0;

      for (i=0; i<=k-1; i++)

      { j=m/2; is=2*is+(m-2*j); m=j;}

      fr[it]=pr[is]; fi[it]=pi[is];

   }

pr[0]=1.0; pi[0]=0.0;

p=6.283185306/(1.0*n);

pr[1]=cos(p); pi[1]=-sin(p);

if (l!=0) pi[1]=-pi[1];

for (i=2; i<=n-1; i++)

   { p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];

      s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);

      pr=p-q; pi=s-p-q;

   }

for (it=0; it<=n-2; it=it+2)

   { vr=fr[it]; vi=fi[it];

      fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];

      fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];

   }

m=n/2; nv=2;

for (l0=k-2; l0>=0; l0--)

   { m=m/2; nv=2*nv;

      for (it=0; it<=(m-1)*nv; it=it+nv)

      for (j=0; j<=(nv/2)-1; j++)

         { p=pr[m*j]*fr[it+j+nv/2];

            q=pi[m*j]*fi[it+j+nv/2];

            s=pr[m*j]+pi[m*j];

            s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);

            poddr=p-q; poddi=s-p-q;

            fr[it+j+nv/2]=fr[it+j]-poddr;

            fi[it+j+nv/2]=fi[it+j]-poddi;

            fr[it+j]=fr[it+j]+poddr;

            fi[it+j]=fi[it+j]+poddi;

         }

   }

if (l!=0)

   for (i=0; i<=n-1; i++)

      { fr=fr/(1.0*n);

      fi=fi/(1.0*n);

      }

if (il!=0)

   for (i=0; i<=n-1; i++)

      { pr=sqrt(fr*fr+fi*fi);

      if (fabs(fr)<0.000001*fabs(fi))

         { if ((fi*fr)>0) pi=90.0;

            else pi=-90.0;

         }

      else

         pi=atan(fi/fr)*360.0/6.283185306;

      }

return;

  }

air23feng · 发表于 2006-4-13 16:47:31

多谢一楼，不知道有没有人在单片机或者ＡＲＭ上实现多

_yuming · 发表于 2006-4-13 17:02:14

这是在PIC18F458里用的FFT程序

/* ************************************************************************

** 函数名:void ADCHD()

** 功能描述:单通道高频采样,数据处理,显示.

*************************************************************************** */

void ADCHD()

{

long D,SHU; //数据.

int n_x,k_x,i; //循环参数.

float Ur,Ui,Urn,Uin; //数据处理中间变量.

ADWORDS=ADCAN+7; //确定MAX197控制字.

for(ADD=0;ADD<64;ADD++) //采样.

{

OUTADC(ADWORDS); //送控制字.

DELAY(0); //延时.

DATA[0][ADD]=INADC(); //读取数据.

}

for(ADD=0;ADD<64;ADD++) //显示波形.

{

D=(int)(40.0*DATA[0][ADD]/4096.0); //取数据.

   LCDPIEX(ADD+64,40-D); //显示.

}

for(n_x=0;n_x<5;n_x++) //计算.

{

Urn=0.0; //实部.

Uin=0.0; //虚部.

for(k_x=0;k_x<32;k_x++) //n_x次谐波.

{

D=DATA[0][k_x]; //取数据计算.

Urn=Urn+D/409.6*cos((2*n_x+1)*(k_x+1)*0.196);

Uin=Uin+D/409.6*sin((2*n_x+1)*(k_x+1)*0.196);

}

Ur=Urn/16.0;                        //

Ui=Uin/16.0; //

SHU=(long)(100*sqrt(Ur*Ur+Ui*Ui));

UI[0][n_x]=SHU; //第n_x次谐波幅值.

UI[0][5]=SHU*SHU+UI[0][5]; //

}

UI[0][5]=(long)sqrt(UI[0][5]); //总幅值.

for(i=0;i<5;i++) //

  {

   SHU=UI[0]*1000; //

     SHU=SHU/(UI[0][5]); //

     UP[0]=SHU; //第i次谐波占有率.

}

SHU=1000*(UI[0][5]-UI[0][0]); //

UP[0][5]=SHU/UI[0][5]; //畸变率.

LCDCLEAN(12,2,126,7); //清除数据显示区.

for(i=0;i<6;i++)

{

OUTNUM(UI[0],1,28,2+i); //显示幅值.

OUTNUM(UP[0],3,56,2+i); //显示占有率.

}

ADWORDS=ADCAN; //还原控制字.

}

air23feng · 发表于 2006-4-13 17:04:36

对_yu-ming 感激不尽！！

_yuming · 发表于 2006-4-13 17:10:55

接3楼

/* ************************************************************************

** 函数名:void XBFX()

** 功能描述:全周波傅立叶积分计算各次谐波的幅值,并返回结果.

*************************************************************************** */

void XBFX()

{

long D,SHU; //数据.

int n_x,k_x,i,n; //循环参数.

float Ur,Ui,Urn,Uin,UIL[2]; //数据处理中间变量.

for(n=0;n<2;n++) //两路信号.

{

for(n_x=0;n_x<5;n_x++) //计算.

{

Urn=0.0; //实部.

Uin=0.0; //虚部.

for(k_x=0;k_x<32;k_x++) //n_x次谐波.

{

D=DATA[CHAN*2+n][k_x]; //取数据计算.

Urn=Urn+D/409.6*cos((2*n_x+1)*(k_x+1)*0.196);

Uin=Uin+D/409.6*sin((2*n_x+1)*(k_x+1)*0.196);

}

Ur=Urn/16.0;                  //

Ui=Uin/16.0;

SHU=(long)(760*sqrt(Ur*Ur+Ui*Ui)); //

UI[n][n_x]=SHU; //

UI[n][5]=SHU*SHU+UI[n][5]; //第n_x次谐波幅值.

if(n_x==0)

UIL[n]=atan(Ui/Ur);    //相位.

}

UI[n][5]=(long)sqrt(UI[n][5]); //总幅值.

for(i=0;i<5;i++) //

    {

      SHU=UI[n]*1000; //

     SHU=SHU/(UI[n][5]); //

     UP[n]=SHU; //第i次谐波占有率.

}

SHU=1000*(UI[n][5]-UI[n][0]); //

UP[n][5]=SHU/UI[n][5]; //畸变率.

}

L[CHAN]=abs((long)(1000*(UIL[0]-UIL[1])-40)); //功率因数角.

S[CHAN]=(long)(UI[0][5]*UI[1][5]); //S.

P[CHAN]=(long)(S[CHAN]*cos(L[CHAN]/1000.0)); //P.

Q[CHAN]=(long)(S[CHAN]*sin(L[CHAN]/1000.0)); //Q.

}

ccdzdpj · 发表于 2006-4-13 17:22:18

不错！顶！

lrzxc · 发表于 2006-4-13 21:36:25

嗯，不错，好帖子。

kinsey · 发表于 2006-4-13 21:55:31

厉害加佩服。

jackrich · 发表于 2006-4-14 08:17:01

正好用到，谢谢！

zhonghua_li · 发表于 2006-8-29 08:23:17

建议不要用浮点数，这样再优化的算法也快不了多少，应采用定点数，我们要的数的范围一般不会那么大的。当然程序写起来也不止于这么简单了。可能要升入到汇编了。

AuToCTRL · 发表于 2006-8-29 08:29:40

To _yu-ming :

你搞工控的？很不错哦。

kfawj · 发表于 2007-3-13 11:01:59

好厉害呀。我要下载下来好好学习一下。

我现在正着急这个傅立叶算法呢！！！

哈哈哈。。。。

y2kloach · 发表于 2007-3-13 11:34:08

虽然暂时用不上,但这么好的算法程序,不顶不行!~~~~~~~~~~~~~~~~~

_yuming · 发表于 2007-3-13 12:42:33

to 10楼，不是所有的地方定点算法都适用，也不是所有的地方浮点算法都适用，要看使用的场合来确定使用哪种算法，很多的地方要求高精度，定点是不能完成的

my_avr · 发表于 2007-3-13 14:59:40

顶一个

sea_19821 · 发表于 2007-3-13 17:52:51

好东西,收藏!

yule · 发表于 2007-3-13 23:02:41

好东西,收藏!

ThinkCell · 发表于 2007-3-14 08:06:19

数学没有学好，还没搞懂计算各次谐波的幅值计算原理。至于如何实现傅立叶积分计算更是不懂，望各位大侠指点一二。谢谢！！

jackrich · 发表于 2007-3-14 08:13:40

支持！

jsnjnzdcyyg · 发表于 2007-3-14 09:11:08

多谢！

starli · 发表于 2007-3-14 13:07:37

可用基2FFT快速算法或基4FFT等,很快的.而且用的存储单元挺少.

我学数字信号处理时学过这个算法.感觉不错.

qsx201 · 发表于 2010-6-1 13:38:26

好东西!

fanwt · 发表于 2010-6-1 13:46:03

不错，目前还没用到

avrwoo · 发表于 2010-6-1 13:49:43

mark

btyang · 发表于 2010-11-23 15:38:53

好东西，收下，谢谢

lan_boy008 · 发表于 2010-11-23 16:19:45

mark

sponge · 发表于 2010-11-23 18:29:36

mark

andrew_dj · 发表于 2010-11-23 19:46:24

z4057 · 发表于 2010-11-23 19:50:26

标记
牛人

ndust · 发表于 2010-11-23 20:26:43

hailin0716 · 发表于 2010-11-23 20:34:56

steven · 发表于 2010-11-23 20:42:55

mark

rlogin · 发表于 2010-11-23 20:43:05

start00 · 发表于 2010-11-23 20:47:23

mark

yqlomg · 发表于 2010-11-23 20:59:14

mark

tu312 · 发表于 2010-11-23 21:03:35

/*
* \file
*
* \brief Header file for FFT library
*
* Copyright (c) 2010 Atmel Corporation. All rights reserved.
*
*
* \page License
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* 3. The name of Atmel may not be used to endorse or promote products derived
* from this software without specific prior written permission.
*
* 4. This software may only be redistributed and used in connection with an
* Atmel AVR product.
*
*  THIS SOFTWARE IS PROVIDED BY ATMEL "AS IS" AND ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT ARE
* EXPRESSLY AND SPECIFICALLY DISCLAIMED. IN NO EVENT SHALL ATMEL BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*/

#include "stdint.h"
#include "stdbool.h"

#define dsp16_t__ int16_t
typedef dsp16_t__ dsp16_t;

//! PI definition also known as the Archimedes' constant
#define PI       (3.141592653589793238462643383279502884197)

/*!
* \brief 32-bit complex signed fixed point type
*/
typedef struct __attribute__ ((__packed__))
{
//! real part
  dsp16_t__ real;
//! imaginary part
  dsp16_t__ imag;
}dsp16_complex_t;

/*! \brief 16-bit fixed point version of the complex FFT algorithm.
*
* \param vout A pointer on a 16-bit complex vector which is the output buffer of this function.
* \param vin A pointer on a 16-bit buffer which is the input buffer of this function.
* \param num It is the size of the input/output vector.
*
*/
extern void FFT (dsp16_complex_t *vout, dsp16_t *vin, uint8_t num);

tu312 · 发表于 2010-11-23 21:06:38

放猛料，Atmel Corporation官方的，这个具体大家懂得
Atmel Corporation官方性能有保证
FFT libraryourdev_599972GCXZTZ.rar(文件大小:3K) (原文件名:libfft.rar)

yywin · 发表于 2010-11-23 21:45:18

学习

yusufu · 发表于 2010-11-23 21:56:11

mark

Adrian · 发表于 2010-11-24 08:44:43

mark

wzyllgx · 发表于 2010-11-24 09:36:38

MARK

xingliu · 发表于 2010-11-24 09:59:40

mark

martial · 发表于 2010-11-24 11:19:46

傅里叶算法，做个标志，好好学习！

10086 · 发表于 2010-11-24 20:22:13

这是好东西，要留

wsdzj · 发表于 2010-12-22 22:30:44

mark

fshunj · 发表于 2010-12-22 22:34:20

好东西

bxzyf · 发表于 2010-12-23 00:13:09

顶顶更健康。

dujun168 · 发表于 2010-12-23 07:51:02

mark

shdjdq · 发表于 2010-12-23 08:40:44

很好的

zengyi703 · 发表于 2010-12-23 08:48:57

mark

liangyurongde · 发表于 2010-12-23 09:20:49

mark

morion · 发表于 2010-12-23 09:23:43

mark了

cunlingwang · 发表于 2010-12-23 10:53:54

学习一下

bjj9217 · 发表于 2010-12-23 10:57:47

mark

jielove2003 · 发表于 2010-12-23 11:18:52

mark

just_interest · 发表于 2010-12-23 12:05:06

ding

deweyled · 发表于 2010-12-23 13:50:59

mark

BINGSHUIHUO · 发表于 2010-12-24 02:52:06

谢谢分享！！

polar · 发表于 2010-12-24 03:03:54

收藏马克一下~

tomhe666 · 发表于 2010-12-24 08:06:57

收藏一下

zhaoghsea · 发表于 2010-12-24 08:23:09

mark

442502587 · 发表于 2010-12-24 08:41:17

收藏马克一下~

jy6715 · 发表于 2010-12-24 08:54:03

mark

jack_yu · 发表于 2010-12-24 09:43:00

好东西.收藏了.谢谢!

liumoz · 发表于 2010-12-24 10:47:12

感激万分。。。。

bbsview · 发表于 2010-12-24 12:52:57

收藏一个

Jerry88 · 发表于 2010-12-24 15:19:14

mark

gc56198 · 发表于 2010-12-24 16:23:35

mark

20061002838 · 发表于 2010-12-24 16:39:08

好东西

wenfeiexe · 发表于 2010-12-24 16:47:48

cool

416446891 · 发表于 2011-1-6 23:06:35

MARK

zzh90513 · 发表于 2011-1-7 01:08:42

学习

jt6245 · 发表于 2011-1-21 11:27:05

mark

zxcao · 发表于 2011-1-21 15:31:14

mark

flyingcys · 发表于 2011-1-24 18:03:54

顶！

wpnx · 发表于 2011-1-24 19:53:45

mark

abcdzhy · 发表于 2011-1-24 23:53:42

mark

onece · 发表于 2011-1-25 08:20:01

mark!mark!

rockgoogle · 发表于 2011-1-25 08:50:59

记号！

lionliu · 发表于 2011-1-25 09:17:59

morion · 发表于 2011-1-25 09:21:05

牛啊收藏了

moen · 发表于 2011-1-25 09:24:13

mark

yanwuxu · 发表于 2011-1-25 09:39:33

快速傅里叶，标记一个^_^

sleet1986 · 发表于 2011-1-25 09:40:15

mark

ypt1980 · 发表于 2011-1-25 10:27:00

/*******************************************************************************
//名称：定点FFT算法程序
//人员：
//日期：
//说明：MEGA128L，16MHZ，IAR420A
//占用周期数：input(5677) Execute(160864) Output(16659)
//16MHZ时ms       0.3548       10.054       1.0412
//fft_r输入信号,范围-255到+255之间，即带符号位9位数据。不在此范围则结果错误。
//输出
*******************************************************************************/
#include "iom128.h"
#include "intrinsics.h"
#include "math.h"

#define  INT8 signed char
#define  INT16  signed int
#define  INT32  signed long

#define  UINT8  unsigned char
#define  UINT16 unsigned int
#define  UINT32 unsigned long
/******************************************************************************/
#define  FFT_N  128 ;             //本程序是128点FFT变换
INT16 fft_r[128]={             //实数采样点，共128点。
//单周期方波
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255
//4周期方波
// 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
//-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
// 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
//-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
// 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
//-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,
// 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
//-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255,-255
//直流恒量
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,
//0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff,0x0ff
};
INT16 fft_temp[128];          //采样点虚部，共128点。
/******************************************************************************/
__flash  UINT16  fft_window[128]={       //汉宁窗
2621,  2639,  2693,  2784,  2910,  3073,  3270,  3502,
3768,  4068,  4401,  4765,  5161,  5587,  6042,  6525,
7036,  7571,  8132,  8715,  9320,  9945,  10588, 11249,
11926, 12616, 13318, 14031, 14753, 15482, 16216, 16954,
17694, 18433, 19171, 19905, 20634, 21356, 22069, 22772,
23462, 24138, 24799, 25443, 26068, 26673, 27256, 27816,
28352, 28862, 29345, 29800, 30226, 30622, 30987, 31319,
31619, 31885, 32117, 32315, 32477, 32603, 32694, 32748,
32766, 32748, 32694, 32603, 32477, 32315, 32117, 31885,
31619, 31319, 30987, 30622, 30226, 29800, 29345, 28862,
28352, 27816, 27256, 26673, 26068, 25443, 24799, 24138,
23462, 22772, 22069, 21356, 20634, 19905, 19171, 18433,
17694, 16954, 16216, 15482, 14753, 14031, 13318, 12616,
11926, 11249, 10588, 9945,  9320,  8715,  8132,  7571,
7036,  6526,  6042,  5587,  5161,  4765,  4401,  4068,
3768,  3502,  3270,  3073,  2910,  2784,  2693,  2639
};
__flash  UINT8 fft_reverse[128]={    //倒位序表
0, 64,32,96, 16,80,48,112,8, 72,40,104,24,88,56,120,
4, 68,36,100,20,84,52,116,12,76,44,108,28,92,60,124,
2, 66,34,98, 18,82,50,114,10,74,42,106,26,90,58,122,
6, 70,38,102,22,86,54,118,14,78,46,110,30,94,62,126,
1, 65,33,97, 17,81,49,113,9, 73,41,105,25,89,57,121,
5, 69,37,101,21,85,53,117,13,77,45,109,29,93,61,125,
3, 67,35,99, 19,83,51,115,11,75,43,107,27,91,59,123,
7, 71,39,103,23,87,55,119,15,79,47,111,31,95,63,127
};
__flash  INT16 fft_wrcos[64]={       //旋转因子实部
32767, 32728, 32609, 32412, 32137, 31785, 31356, 30852,
30273, 29621, 28898, 28105, 27245, 26319, 25329, 24279,
23170, 22005, 20787, 19519, 18204, 16846, 15446, 14010,
12539, 11039, 9512,  7962,  6393,  4808,  3212,  1608,
0,    -1608, -3212, -4808, -6393, -7962, -9512, -11039,
-12539,-14010,-15446,-16846,-18204,-19519,-20787,-22005,
-23170,-24279,-25329,-26319,-27245,-28105,-28898,-29621,
-30273,-30852,-31356,-31785,-32137,-32412,-32609,-32728
};
__flash  INT16  ttt_wisin[64]={       //旋转因子虚部
0,    1608,  3212,  4808,  6393,  7962,  9512,  11039,
12539, 14010, 15446, 16846, 18204, 19519, 20787, 22005,
23170, 24279, 25329, 26319, 27245, 28105, 28898, 29621,
30273, 30852, 31356, 31785, 32137, 32412, 32609, 32728,
32767, 32728, 32609, 32412, 32137, 31785, 31356, 30852,
30273, 29621, 28898, 28105, 27245, 26319, 25329, 24279,
23170, 22005, 20787, 19519, 18204, 16846, 15446, 14010,
12539, 11039, 9512,  7962,  6393,  4808,  3212,  1608
};
/******************************************************************************/
void fft_real( void );                   //fft变换函数
UINT16  sqrt_16( UINT32 M );                //开平方函数
extern  UINT32  fft_muls16( INT16 mula, INT16 mulb ); //16x16位乘法运算函数

volatile UINT16 f[42],ff;                   //存放结果
/******************************************************************************/
void main( void )
{

  while(1){
fft_real( );                         //fft变换函数

while(1);                            //变换后结果取代原来数据，所以仅能运行1次
  }
}
/******************************************************************************/
//名称： fft_real（）定点数快速傅立叶变换函数
//参数：无入口和出口参数
//说明：计算42次谐波
/******************************************************************************/
void fft_real( void )
{
  UINT8 i,j,k;
  UINT8 l,b,p;
  UINT8 m;
  INT16 tr,ti;
  INT16 rcos_isin,rsin_icos;
  INT16 *pa,*pc;

  pa= fft_r;
  pc= fft_temp;
  for(i=0; i<128; i++)
   {
   *pc =  *pa;  //不加窗
   //*pc = ((INT32)(*pa)*(fft_window))>>15;  //((UINT32)(*pa)*(*pb))>>15; //加窗
   pa++;  pc++;
   }
  pa= fft_r;
  pc= fft_temp;
  for(i=0; i<128; i++)
   {
   pa = pc[fft_reverse]; //倒位序
   pc[fft_reverse] = 0;       //虚部为0
   }
  b= 1;                            //b=2^l =1 ,2 ,4 ,8,16,32,64
  i= 6;                            //p=2^(6-l)=64(6),32(5),16(4),8(3),4(2),2(1),1(0)
  for(l=0; l<7; l++)                //级数
   {
   for(j=0; j<b; j++)          //旋转因子指数
      {
      p= (j<<i);                //p=2^(6-l)=64(6),32(5),16(4),8(3),4(2),2(1),1(0)
      for(k=j; k<128; )       //蝶形
         {
            m = k+b;
            tr=fft_r[k];
            ti=fft_temp[k];
                                                //表中值最大为32768, 32768/(2^15)=1.
            rcos_isin= ((INT32)fft_r[m]* fft_wrcos[p]+(INT32)fft_temp[m]* ttt_wisin[p])>>15;
            rsin_icos= ((INT32)fft_r[m]* ttt_wisin[p]-(INT32)fft_temp[m]* fft_wrcos[p])>>15;
            //rcos_isin= (fft_muls16(fft_r[m], fft_wrcos[p]) + fft_muls16(fft_temp[m], ttt_wisin[p]))>>15;
            //rsin_icos= (fft_muls16(fft_r[m], ttt_wisin[p]) + fft_muls16(fft_temp[m], fft_wrcos[p]))>>15;
            fft_r[k]    =(tr+ rcos_isin);    //为了保证TR，TI不溢出，输入信号位数9位。
            fft_temp[k]  =(ti- rsin_icos);
            fft_r[m]    =(tr- rcos_isin);
            fft_temp[m]  =(ti+ rsin_icos);

            k = m+b;    //k=k+2*b
         }             //for k
      }                //for j
   i--;                //p=2^(6-l)=64,32,16,8,4 ,2 ,1
   b<<=1;             //b=2^l =1 ,2 ,4 ,8,16,32,64
   }                   //for l
for(i=0; i<42; i++)       //求得的谐波频率幂次小于点数的1/3。
{
   f= sqrt_16((INT32)((INT32)fft_r*fft_r+ (INT32)fft_temp*fft_temp) );
   // f= sqrt_16( fft_muls16(fft_r, fft_r) + fft_muls16(fft_temp, fft_temp) );
}
  ff=f[0];  //f[4]=20774; f[12]=7010; f[20]=4318; f[28]=3208; f[36]=2635;
  //直流量基波       3次       5次       7次       9次
}
/******************************************************************************/
//名称：sqrt_16（）
//参数：M入口参数32位，N出口参数16位
//说明：开平方函数
/******************************************************************************/
UINT16 sqrt_16( UINT32 M )
{
UINT16  N, i;
UINT32  tmp, ttp;       // 结果、循环计数
if(M == 0)             // 被开方数，开方结果也为0
   return 0;
N = 0;
tmp = (M >> 30);       // 获取最高位：B[m-1]
M <<= 2;
if(tmp > 1)             // 最高位为1
{
   N ++;             // 结果当前位为1，否则为默认的0
   tmp -= N;
}
for(i=15; i>0; i--)    // 求剩余的15位
{
   N <<= 1;          // 左移一位
   tmp <<= 2;
   tmp += (M >> 30); // 假设
   ttp = N;
   ttp = (ttp<<1)+1;
   M <<= 2;
   if (tmp >= ttp)    // 假设成立
      {
      tmp -= ttp;
      N ++;
      }
}
return N;
}
/******************************************************************************/

almasy · 发表于 2011-1-25 10:42:51

mark

geniusjia · 发表于 2011-1-25 14:31:06

mark

gxy508 · 发表于 2011-1-25 15:28:21

mark

hecb999 · 发表于 2011-1-25 16:44:48

FFT mark

fesu02 · 发表于 2011-1-25 18:44:42

zgcumt · 发表于 2011-1-26 08:57:21

再次mark

myhonour · 发表于 2011-1-26 11:50:03

MARK

zhangjinxing · 发表于 2011-1-26 12:07:04

mark

xiaomage_2000 · 发表于 2011-1-26 12:37:45

我一直不明白这个，打个比方我一个8位的AD读入比如音频电压。是如何把这个8位数据运算成频谱显示的？

bingshuihuo888 · 发表于 2011-1-26 15:25:43

谢谢

hizzx · 发表于 2011-1-26 16:33:15

做个记号。

qsx201 · 发表于 2011-1-26 16:50:47

mark

ringan865 · 发表于 2011-1-26 19:02:07

mark

elecom · 发表于 2011-1-26 20:13:38

mark

寻求快速傅立叶算法的Ｃ语言实现

阿莫论坛20周年了！感谢大家的支持与爱护！！