寻求快速傅立叶算法的Ｃ语言实现 (amobbs.com 阿莫电子技术论坛) -

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=pr; fi=pi;

   }

pr=1.0; pi=0.0;

p=6.283185306/(1.0*n);

pr=cos(p); pi=-sin(p);

if (l!=0) pi=-pi;

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

   { p=pr*pr; q=pi*pi;

   s=(pr+pi)*(pr+pi);

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

   }

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

   { vr=fr; vi=fi;

   fr=vr+fr; fi=vi+fi;

   fr=vr-fr; fi=vi-fi;

   }

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*fr;

         q=pi*fi;

         s=pr+pi;

         s=s*(fr+fi);

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

         fr=fr-poddr;

         fi=fi-poddi;

         fr=fr+poddr;

         fi=fi+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=INADC(); //读取数据.

}

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

{

D=(int)(40.0*DATA/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; //取数据计算.

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=SHU; //第n_x次谐波幅值.

UI=SHU*SHU+UI; //

}

UI=(long)sqrt(UI); //总幅值.

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

{

   SHU=UI*1000; //

   SHU=SHU/(UI); //

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

}

SHU=1000*(UI-UI); //

UP=SHU/UI; //畸变率.

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

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

{

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

OUTNUM(UP,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; //数据处理中间变量.

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; //取数据计算.

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=SHU; //

UI=SHU*SHU+UI; //第n_x次谐波幅值.

if(n_x==0)

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

}

UI=(long)sqrt(UI); //总幅值.

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

{

      SHU=UI*1000; //

   SHU=SHU/(UI); //

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

}

SHU=1000*(UI-UI); //

UP=SHU/UI; //畸变率.

}

L=abs((long)(1000*(UIL-UIL)-40)); //功率因数角.

S=(long)(UI*UI); //S.

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

Q=(long)(S*sin(L/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

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

./emotion/em026.gif

_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"

#defineINT8 signed char
#defineINT16signed int
#defineINT32signed long

#defineUINT8unsigned char
#defineUINT16 unsigned int
#defineUINT32 unsigned long
/******************************************************************************/
#defineFFT_N128 ;             //本程序是128点FFT变换
INT16 fft_r={             //实数采样点，共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点。
/******************************************************************************/
__flashUINT16fft_window={       //汉宁窗
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
};
__flashUINT8 fft_reverse={    //倒位序表
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
};
__flashINT16 fft_wrcos={       //旋转因子实部
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
};
__flashINT16ttt_wisin={       //旋转因子虚部
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变换函数
UINT16sqrt_16( UINT32 M );                //开平方函数
externUINT32fft_muls16( INT16 mula, INT16 mulb ); //16x16位乘法运算函数

volatile UINT16 f,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]; //倒位序
   pc] = 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;
         ti=fft_temp;
                                             //表中值最大为32768, 32768/(2^15)=1.
         rcos_isin= ((INT32)fft_r* fft_wrcos+(INT32)fft_temp* ttt_wisin)>>15;
         rsin_icos= ((INT32)fft_r* ttt_wisin-(INT32)fft_temp* fft_wrcos)>>15;
         //rcos_isin= (fft_muls16(fft_r, fft_wrcos) + fft_muls16(fft_temp, ttt_wisin))>>15;
         //rsin_icos= (fft_muls16(fft_r, ttt_wisin) + fft_muls16(fft_temp, fft_wrcos))>>15;
         fft_r =(tr+ rcos_isin);    //为了保证TR，TI不溢出，输入信号位数9位。
         fft_temp=(ti- rsin_icos);
         fft_r =(tr- rcos_isin);
         fft_temp=(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;//f=20774; f=7010; f=4318; f=3208; f=2635;
//直流量基波       3次       5次       7次       9次
}
/******************************************************************************/
//名称：sqrt_16（）
//参数：M入口参数32位，N出口参数16位
//说明：开平方函数
/******************************************************************************/
UINT16 sqrt_16( UINT32 M )
{
UINT16N, i;
UINT32tmp, ttp;    // 结果、循环计数
if(M == 0)             // 被开方数，开方结果也为0
return 0;
N = 0;
tmp = (M >> 30);       // 获取最高位：B
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

FFTmark

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

页: [1] 2

amobbs.com 阿莫电子技术论坛's Archiver

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