回复: 183
寻求快速傅立叶算法的C语言实现
(476465702)
出0入0汤圆
电梯直达
发表于 2006-4-13 16:29:22
|
只看该作者
|倒序浏览
|阅读模式
(476465270)
出0入0汤圆
发表于 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;
}
(476464613)
出0入0汤圆
楼主|
发表于 2006-4-13 16:47:31
|
只看该作者
多谢一楼,不知道有没有人在单片机或者ARM上实现多
(476463730)
出0入0汤圆
发表于 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; //还原控制字.
}
(476463588)
出0入0汤圆
楼主|
发表于 2006-4-13 17:04:36
|
只看该作者
对_yu-ming 感激不尽!!
(476463209)
出0入0汤圆
发表于 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.
}
(476462526)
出0入0汤圆
发表于 2006-4-13 17:22:18
|
只看该作者
不错!顶!
(476447279)
出0入0汤圆
发表于 2006-4-13 21:36:25
|
只看该作者
嗯,不错,好帖子。
(476446133)
出0入0汤圆
发表于 2006-4-13 21:55:31
|
只看该作者
厉害加佩服。
(476408843)
出0入0汤圆
发表于 2006-4-14 08:17:01
|
只看该作者
正好用到,谢谢!
发表于 2006-8-29 08:23:17
|
只看该作者
建议不要用浮点数,这样再优化的算法也快不了多少,应采用定点数,我们要的数的范围一般不会那么大的。当然程序写起来也不止于这么简单了。可能要升入到汇编了。
(464571284)
出0入0汤圆
发表于 2006-8-29 08:29:40
|
只看该作者
To _yu-ming :
你搞工控的?很不错哦。
积分
精华汤圆新注册者
{*}
(447627745)
出0入0汤圆
发表于 2007-3-13 11:01:59
|
只看该作者
好厉害呀。我要下载下来好好学习一下。
我现在正着急这个傅立叶算法呢!!!
哈哈哈。。。。
(447625816)
出0入0汤圆
发表于 2007-3-13 11:34:08
|
只看该作者
虽然暂时用不上,但这么好的算法程序,不顶不行!~~~~~~~~~~~~~~~~~
(447621711)
出0入0汤圆
发表于 2007-3-13 12:42:33
|
只看该作者
to 10楼,不是所有的地方定点算法都适用,也不是所有的地方浮点算法都适用,要看使用的场合来确定使用哪种算法,很多的地方要求高精度,定点是不能完成的
(447613484)
出0入0汤圆
发表于 2007-3-13 14:59:40
|
只看该作者
顶一个
(447603093)
出0入0汤圆
发表于 2007-3-13 17:52:51
|
只看该作者
好东西,收藏!
当前离线
精华汤圆游客
{*}
(447584503)
出0入0汤圆
发表于 2007-3-13 23:02:41
|
只看该作者
好东西,收藏!
(447551885)
出0入0汤圆
发表于 2007-3-14 08:06:19
|
只看该作者
数学没有学好,还没搞懂计算各次谐波的幅值计算原理。至于如何实现傅立叶积分计算更是不懂,望各位大侠指点一二。谢谢!!
(447551444)
出0入0汤圆
发表于 2007-3-14 08:13:40
|
只看该作者
支持!
(447547996)
出0入0汤圆
发表于 2007-3-14 09:11:08
|
只看该作者
多谢!
(447533807)
出0入0汤圆
发表于 2007-3-14 13:07:37
|
只看该作者
可用基2FFT快速算法或基4FFT等,很快的.而且用的存储单元挺少.
我学数字信号处理时学过这个算法.感觉不错.
积分
精华汤圆新手上路
{*}
(346011958)
出0入0汤圆
发表于 2010-6-1 13:38:26
|
只看该作者
好东西!
(346011501)
出0入0汤圆
发表于 2010-6-1 13:46:03
|
只看该作者
不错,目前还没用到
(346011281)
出0入0汤圆
发表于 2010-6-1 13:49:43
|
只看该作者
mark
积分
精华汤圆注册会员
{*}
(330884731)
出0入0汤圆
发表于 2010-11-23 15:38:53
|
只看该作者
好东西,收下,谢谢
(330882279)
出0入0汤圆
发表于 2010-11-23 16:19:45
|
只看该作者
mark
积分
精华汤圆注册会员
{*}
(330874488)
出0入0汤圆
发表于 2010-11-23 18:29:36
|
只看该作者
mark
(330869880)
出0入0汤圆
发表于 2010-11-23 19:46:24
|
只看该作者
收藏
精华汤圆中级会员
{*}
(330869638)
出0入0汤圆
发表于 2010-11-23 19:50:26
|
只看该作者
标记
牛人
(330867461)
出0入0汤圆
发表于 2010-11-23 20:26:43
|
只看该作者
jh
(330866968)
出0入0汤圆
发表于 2010-11-23 20:34:56
|
只看该作者
顶
(330866489)
出0入0汤圆
发表于 2010-11-23 20:42:55
|
只看该作者
mark
(330866479)
出0入0汤圆
发表于 2010-11-23 20:43:05
|
只看该作者
~
(330866221)
出0入0汤圆
发表于 2010-11-23 20:47:23
|
只看该作者
mark
(330865510)
出0入0汤圆
发表于 2010-11-23 20:59:14
|
只看该作者
mark
积分
精华汤圆注册会员
{*}
(330865249)
出0入0汤圆
发表于 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);
积分
精华汤圆注册会员
{*}
(330865066)
出0入0汤圆
发表于 2010-11-23 21:06:38
|
只看该作者
放猛料,Atmel Corporation官方的,这个具体大家懂得
Atmel Corporation官方性能有保证
FFT libraryourdev_599972GCXZTZ.rar(文件大小:3K) (原文件名:libfft.rar)
精华汤圆中级会员
{*}
(330862746)
出0入0汤圆
发表于 2010-11-23 21:45:18
|
只看该作者
学习
精华汤圆VIP
{*}
(330862093)
出0入0汤圆
发表于 2010-11-23 21:56:11
|
只看该作者
mark
(330823181)
出0入0汤圆
发表于 2010-11-24 08:44:43
|
只看该作者
mark
(330820066)
出0入0汤圆
发表于 2010-11-24 09:36:38
|
只看该作者
MARK
(330818684)
出0入0汤圆
发表于 2010-11-24 09:59:40
|
只看该作者
mark
积分
精华汤圆新注册者
{*}
(330813878)
出0入0汤圆
发表于 2010-11-24 11:19:46
|
只看该作者
傅里叶算法,做个标志,好好学习!
精华汤圆中级会员
{*}
(330781331)
出0入0汤圆
发表于 2010-11-24 20:22:13
|
只看该作者
这是好东西,要留
积分
精华汤圆新手上路
{*}
(328354420)
出0入0汤圆
发表于 2010-12-22 22:30:44
|
只看该作者
mark
(328354204)
出0入0汤圆
发表于 2010-12-22 22:34:20
|
只看该作者
好东西
精华汤圆高级会员
{*}
(328348275)
出0入0汤圆
发表于 2010-12-23 00:13:09
|
只看该作者
顶顶更健康。
(328320802)
出0入0汤圆
发表于 2010-12-23 07:51:02
|
只看该作者
mark
(328317820)
出0入0汤圆
发表于 2010-12-23 08:40:44
|
只看该作者
很好的
(328317327)
出0入0汤圆
发表于 2010-12-23 08:48:57
|
只看该作者
mark
发表于 2010-12-23 09:20:49
|
只看该作者
mark
(328315241)
出0入0汤圆
发表于 2010-12-23 09:23:43
|
只看该作者
mark了
(328309830)
出0入0汤圆
发表于 2010-12-23 10:53:54
|
只看该作者
学习一下
(328309597)
出0入0汤圆
发表于 2010-12-23 10:57:47
|
只看该作者
mark
发表于 2010-12-23 11:18:52
|
只看该作者
mark
发表于 2010-12-23 12:05:06
|
只看该作者
ding
(328299205)
出0入0汤圆
发表于 2010-12-23 13:50:59
|
只看该作者
mark
(328252338)
出0入0汤圆
发表于 2010-12-24 02:52:06
|
只看该作者
谢谢分享!!
精华汤圆注册会员
{*}
(328251630)
出0入0汤圆
发表于 2010-12-24 03:03:54
|
只看该作者
收藏 马克一下~
(328233447)
出0入0汤圆
发表于 2010-12-24 08:06:57
|
只看该作者
收藏一下
(328232475)
出0入0汤圆
发表于 2010-12-24 08:23:09
|
只看该作者
mark
(328231387)
出0入12汤圆
发表于 2010-12-24 08:41:17
|
只看该作者
收藏 马克一下~
(328230621)
出0入0汤圆
发表于 2010-12-24 08:54:03
|
只看该作者
mark
(328227684)
出0入0汤圆
发表于 2010-12-24 09:43:00
|
只看该作者
好东西.收藏了.谢谢!
积分
精华汤圆新手上路
{*}
(328223832)
出0入0汤圆
发表于 2010-12-24 10:47:12
|
只看该作者
感激万分。。。。
(328216287)
出0入0汤圆
发表于 2010-12-24 12:52:57
|
只看该作者
收藏一个
积分
精华汤圆注册会员
{*}
(328207510)
出0入0汤圆
发表于 2010-12-24 15:19:14
|
只看该作者
mark
(328203649)
出0入0汤圆
发表于 2010-12-24 16:23:35
|
只看该作者
mark
发表于 2010-12-24 16:39:08
|
只看该作者
好东西
(328202196)
出0入0汤圆
发表于 2010-12-24 16:47:48
|
只看该作者
cool
(327056269)
出0入0汤圆
发表于 2011-1-6 23:06:35
|
只看该作者
MARK
(327048942)
出0入0汤圆
发表于 2011-1-7 01:08:42
|
只看该作者
学习
(325802239)
出0入0汤圆
发表于 2011-1-21 11:27:05
|
只看该作者
mark
精华汤圆中级会员
{*}
(325787590)
出0入0汤圆
发表于 2011-1-21 15:31:14
|
只看该作者
mark
(325519230)
出0入0汤圆
发表于 2011-1-24 18:03:54
|
只看该作者
顶!
当前离线
精华汤圆中级会员
{*}
(325512639)
出0入0汤圆
发表于 2011-1-24 19:53:45
|
只看该作者
mark
(325498242)
出0入0汤圆
发表于 2011-1-24 23:53:42
|
只看该作者
mark
(325467863)
出0入0汤圆
发表于 2011-1-25 08:20:01
|
只看该作者
mark!mark!
(325466005)
出0入0汤圆
发表于 2011-1-25 08:50:59
|
只看该作者
记号!
(325464385)
出0入0汤圆
发表于 2011-1-25 09:17:59
|
只看该作者
m
(325464199)
出0入0汤圆
发表于 2011-1-25 09:21:05
|
只看该作者
牛啊 收藏了
当前离线
精华汤圆VIP++
{*}
(325464011)
出0入0汤圆
发表于 2011-1-25 09:24:13
|
只看该作者
mark
(325463091)
出0入0汤圆
发表于 2011-1-25 09:39:33
|
只看该作者
快速傅里叶,标记一个^_^
(325463049)
出0入0汤圆
发表于 2011-1-25 09:40:15
|
只看该作者
mark
积分
精华汤圆新手上路
{*}
(325460244)
出0入0汤圆
发表于 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
{
p= (j<
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;
}
/******************************************************************************/
(325459293)
出0入0汤圆
发表于 2011-1-25 10:42:51
|
只看该作者
mark
(325445598)
出0入0汤圆
发表于 2011-1-25 14:31:06
|
只看该作者
mark
(325442163)
出0入0汤圆
发表于 2011-1-25 15:28:21
|
只看该作者
mark
(325437576)
出0入0汤圆
发表于 2011-1-25 16:44:48
|
只看该作者
FFT mark
积分
精华汤圆新手上路
{*}
(325430382)
出0入0汤圆
发表于 2011-1-25 18:44:42
|
只看该作者
M
积分
精华汤圆注册会员
{*}
(325379223)
出0入0汤圆
发表于 2011-1-26 08:57:21
|
只看该作者
再次mark
(325368861)
出0入0汤圆
发表于 2011-1-26 11:50:03
|
只看该作者
MARK
发表于 2011-1-26 12:07:04
|
只看该作者
mark
发表于 2011-1-26 12:37:45
|
只看该作者
我一直不明白这个,打个比方我一个8位的AD读入比如音频电压。是如何把这个8位数据运算成频谱显示的?
发表于 2011-1-26 15:25:43
|
只看该作者
谢谢
积分
精华汤圆注册会员
{*}
(325351869)
出0入0汤圆
发表于 2011-1-26 16:33:15
|
只看该作者
做个记号。
积分
精华汤圆新手上路
{*}
(325350817)
出0入0汤圆
发表于 2011-1-26 16:50:47
|
只看该作者
mark
(325342937)
出0入0汤圆
发表于 2011-1-26 19:02:07
|
只看该作者
mark
(325338646)
出0入0汤圆
发表于 2011-1-26 20:13:38
|
只看该作者
mark