寻求快速傅立叶算法的C语言实现
如题,望大虾点拨! #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;
} 多谢一楼,不知道有没有人在单片机或者ARM上实现多 这是在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; //还原控制字.
} 对_yu-ming 感激不尽!! 接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.
} 不错!顶! 嗯,不错,好帖子。 厉害加佩服。 正好用到,谢谢! 建议不要用浮点数,这样再优化的算法也快不了多少,应采用定点数,我们要的数的范围一般不会那么大的。当然程序写起来也不止于这么简单了。可能要升入到汇编了。 To _yu-ming :
你搞工控的?很不错哦。 好厉害呀。我要下载下来好好学习一下。
我现在正着急这个傅立叶算法呢!!!
哈哈哈。。。。 虽然暂时用不上,但这么好的算法程序,不顶不行!~~~~~~~~~~~~~~~~~
./emotion/em026.gif to 10楼,不是所有的地方定点算法都适用,也不是所有的地方浮点算法都适用,要看使用的场合来确定使用哪种算法,很多的地方要求高精度,定点是不能完成的 顶一个 好东西,收藏! 好东西,收藏! 数学没有学好,还没搞懂计算各次谐波的幅值计算原理。至于如何实现傅立叶积分计算更是不懂,望各位大侠指点一二。谢谢!! 支持! 多谢! 可用基2FFT快速算法或基4FFT等,很快的.而且用的存储单元挺少.
我学数字信号处理时学过这个算法.感觉不错. 好东西! 不错,目前还没用到 mark 好东西,收下,谢谢 mark mark 收藏 标记
牛人 jh 顶 mark ~ mark mark /*
* \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); 放猛料,Atmel Corporation官方的,这个具体大家懂得
Atmel Corporation官方性能有保证
FFT libraryourdev_599972GCXZTZ.rar(文件大小:3K) (原文件名:libfft.rar) 学习 mark mark MARK mark 傅里叶算法,做个标志,好好学习! 这是好东西,要留 mark 好东西 顶顶更健康。 mark 很好的 mark mark mark了 学习一下 mark mark ding mark 谢谢分享!! 收藏 马克一下~ 收藏一下 mark 收藏 马克一下~ mark 好东西.收藏了.谢谢! 感激万分。。。。 收藏一个 mark mark 好东西 cool MARK 学习 mark mark 顶! mark mark mark!mark! 记号! m 牛啊 收藏了 mark 快速傅里叶,标记一个^_^ mark /*******************************************************************************
//名称:定点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;
}
/******************************************************************************/ mark mark mark FFTmark M 再次mark MARK mark 我一直不明白这个,打个比方我一个8位的AD读入比如音频电压。是如何把这个8位数据运算成频谱显示的? 谢谢 做个记号。 mark mark mark
页:
[1]
2