STM32音乐频谱分析 快速傅里叶变换FFT资料与源码分享

内有FFT算法资料，程式及程式详细说明等



全部资料51hei下载地址：
音乐频谱分析fft.rar (19.68 MB, 下载次数: 555)

2018-10-10 11:42 上传

点击文件名下载附件

FFT位翻转程序效率：
在单片机控制程序中，往往会用到位翻转程序，例如点阵的控制，图形的处理，FFT运算等。那么，在C语言中如何才能写出高效率的程序呢？今日在keil的论坛中看到有网友提及这个程序，又在论坛搜索了一下，将老外写的，网友写的，我自己写的程序做了一个全方位的测试，结果如下所示：
首先是老外的程序：
作者：Concepcion Marco Valero
#include <reg52.h>
unsigned char mr;
unsigned char invertir_byte (mr) {
mr = (mr & 0x0F) << 4 | (mr & 0xF0) >> 4;
mr = (mr & 0x33) << 2 | (mr & 0xCC) >> 2;
mr = (mr & 0x55) << 1 | (mr & 0xAA) >> 1;
return (mr);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
Program Size: data=10.0 xdata=0 code=123
完成位交换需要 121 个时钟周期。
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第二个程序：我写的
#include <reg52.h>
unsigned char mr;
unsigned char invertir_byte (mr) {
unsigned char temp;
if(mr&0x80){temp=temp|0x01;}
if(mr&0x40){temp=temp|0x02;}
if(mr&0x20){temp=temp|0x04;}
if(mr&0x10){temp=temp|0x08;}
if(mr&0x08){temp=temp|0x10;}
if(mr&0x04){temp=temp|0x20;}
if(mr&0x02){temp=temp|0x40;}
if(mr&0x01){temp=temp|0x80;}
return (temp);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
Program Size: data=10.0 xdata=0 code=85
完成位交换需要 42 个时钟周期。
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■

第三个程序：还是我写的
#include <reg52.h>
unsigned char mr;
unsigned char invertir_byte (mr) {
bit tempb;
unsigned char count,temp;
for(count=8;count;count--)
{
tempb=mr&0x01;
mr>>=1;
temp<<=1;
temp=temp|tempb;
}
return (temp);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
Program Size: data=12.1 xdata=0 code=64
完成位交换需要 175 个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第三个程序：还是我写的
#include <reg52.h>
unsigned char mr;
unsigned char invertir_byte (mr) {
bit tempb;
unsigned char count,temp;
for(count=8;count;count--)
{
tempb=mr&0x01;
mr>>=1;
temp<<=1;
temp=temp|tempb;
}
return (temp);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
Program Size: data=12.1 xdata=0 code=64
完成位交换需要 175 个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第四个程序：还是我写的
#include <reg52.h>
unsigned char bdata temp;
sbit D0=temp^0;
sbit D1=temp^1;
sbit D2=temp^2;
sbit D3=temp^3;
sbit D4=temp^4;
sbit D5=temp^5;
sbit D6=temp^6;
sbit D7=temp^7;
unsigned char invertir_byte (unsigned char mr)
{
D7=mr&0x01;
D6=mr&0x02;
D5=mr&0x04;
D4=mr&0x08;
D3=mr&0x10;
D2=mr&0x20;
D1=mr&0x40;
D0=mr&0x80;
return (temp);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
Program Size: data=10.0 xdata=0 code=59
完成位交换需要 35个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第五个程序：Jon Ward
##include <reg52.h>
unsigned char bdata src;
sbit S0=src^0;
sbit S1=src^1;
sbit S2=src^2;
sbit S3=src^3;
sbit S4=src^4;
sbit S5=src^5;
sbit S6=src^6;
sbit S7=src^7;
unsigned char bdata dst;
sbit D0=dst^0;
sbit D1=dst^1;
sbit D2=dst^2;
sbit D3=dst^3;
sbit D4=dst^4;
sbit D5=dst^5;
sbit D6=dst^6;
sbit D7=dst^7;
unsigned char invertir_byte (unsigned char mr)
{
src=mr;
D0=S7;
D1=S6;
D2=S5;
D3=S4;
D4=S3;
D5=S2;
D6=S1;
D7=S0;
return(dst);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
//cost 35 machine cycle
//Program Size: data=11.0 xdata=0 code=61
完成位交换需要 35个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第六个程序：来自Ourdev论坛的网友
#include <reg52.h>
unsigned char invertir_byte (unsigned char val)
{
unsigned char dat_b ,i;
dat_b=0x00;
for(i=0;i<=7;i++)
{
dat_b=dat_b|((val>>i)&0x01);
if(i==7)break;
dat_b=dat_b<<1;
}
val=dat_b;
return(val);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}

287 cycle
Program Size: data=9.0 xdata=0 code=57
完成位交换需要 287个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第七个程序：来自ourdev论坛的网友
#include <reg52.h>
unsigned char code tab[16]={0x00,0x08,0x04,0x0c,0x02,0x0a,0x06,0x0e,
0x01,0x09,0x05,0x0d,0x03,0x0b,0x07,0x0f};
unsigned char invertir_byte (unsigned char dat)
{
dat = tab[(dat & 0xf0)>>4] | (tab[dat & 0x0f]<<4);
return dat;
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}
//cost 26 machine cycle
//Program Size: data=9.0 xdata=0 code=63
完成位交换需要 26 个时钟周期
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
第八个程序：来自ourdev网友
#include<AT89X51.H>
unsigned char byte_bit_swap(unsigned char a)
{
a = ((a & 0x0F) << 4) | ((a & 0xF0) >> 4);
a = ((a << 2) & 0xcc) | ((a>> 2) & 0x33);
a = ((a << 1) & 0xaa) | ((a>> 1) & 0x55);
return(a);
}
void main(void)
{
while(1)
{
P1=byte_bit_swap(0x33);
}
}
Program Size: data=9.0 xdata=0 code=66
完成位交换需要 37 个时钟周期
　



// 快速傅里 叶变换FFT的C语言算法彻底研究
         // LED音乐频谱显示的核心算法就是快速傅里叶变换，FFT的理解和编程还是比较难的，特地撰写此文分享一下研究成果。
// 一、彻底理解傅里叶变换   
// 快速傅里叶变换（Fast Fourier Transform)是离散傅里叶变换的一种快速算法，简称FFT，通过FFT可以将一个信号从时域变换到频域。
// 模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。
// 假设采样频率为fs，采样点数为N，那么FFT结果就是一个N点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/N。
// 举例说明：用1kHz的采样频率采样128点，则FFT结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 Hz。
// 这个频率点的幅值为：该点复数的模值除以N/2（n=1时是直流分量，其幅值是该点的模值除以N）。
// 二、傅里叶变换的C语言编程
// 1、对于快速傅里叶变换FFT，第一个要解决的问题就是码位倒序。
    // 假设一个N点的输入序列，那么它的序号二进制数位数就是t=log2N.
    // 码位倒序要解决两个问题：①将t位二进制数倒序；②将倒序后的两个存储单元进行交换。

// ①将t=3位二进制数倒序

// ②将倒序后的两个存储单元进行交换
    // 如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢？利用C语言的移位功能！
// 程序如下，我不多说，看不懂者智商一定在180以下！
// 复数类型定义及其运算
#define N 64      //64点
#define log2N 6   //log2N=6
/*复数类型*/
typedef struct
{
float real;
float img;
}complex;
complex xdata x[N]; //输入序列
/*复数加法*/
complex add(complex a,complex b)
{
complex c;
c.real=a.real+b.real;
c.img=a.img+b.img;
return c;
}
/*复数减法*/
complex sub(complex a,complex b)
{
complex c;
c.real=a.real-b.real;
c.img=a.img-b.img;
return c;
}
/*复数乘法*/
complex mul(complex a,complex b)
{
complex c;
c.real=a.real*b.real - a.img*b.img;
c.img=a.real*b.img + a.img*b.real;
return c;
}
/***码位倒序函数***/
void Reverse(void)
{
unsigned int i,j,k;
unsigned int t;
complex temp;//临时交换变量
for(i=0;i<N;i++)//从第0个序号到第N-1个序号
{
  k=i;//当前第i个序号
  j=0;//存储倒序后的序号，先初始化为0
  for(t=0;t<log2N;t++)//共移位t次，其中log2N是事先宏定义算好的
  {
   j<<=1;
   j|=(k&1);//j左移一位然后加上k的最低位
   k>>=1;//k右移一位，次低位变为最低位
  }
  if(j>i)//如果倒序后大于原序数，就将两个存储单元进行交换（判断j>i是为了防止重复交换）
  {
   temp=x[ i];
   x[ i]=x[j];
   x[j]=temp;
  }
}
}
2、第二个要解决的问题就是蝶形运算
  
①第1级（第1列）每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，
第m级蝶形运算，每个蝶形的两节点“距离”L=2m-1。
②对于16点的FFT，第1级有16组蝶形，每组有1个蝶形；第2级有4组蝶形，每组有2个蝶形；第3级有2组蝶形，每组有4个蝶形；第4级有1组蝶形，每组有8个蝶形。由此可推出，
对于N点的FFT，第m级有N/2L组蝶形，每组有L=2m-1个蝶形。
③旋转因子 的确定
以16点FFT为例，第m级第k个旋转因子为 ，其中k=0～2m-1-1，即第m级共有2m-1个旋转因子，根据旋转因子的可约性， ，所以第m级第k个旋转因子为 ，其中k=0～2m-1-1。
为提高FFT的运算速度，我们可以事先建立一个旋转因子数组，然后通过查表法来实现。
complex code WN[N]=//旋转因子数组
{ //为节省CPU计算时间，旋转因子采用查表处理
  //★根据实际FFT的点数N，该表数据需自行修改
  //以下结果通过Excel自动生成
  //  WN[k].real=cos(2*PI/N*k);
  //  WN[k].img=-sin(2*PI/N*k);
  {1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
  {0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
  {0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
  {0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
  {0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
  {-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
  {-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
  {-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
  {-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
  {-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
  {-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
  {-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
  {0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
  {0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
  {0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
  {0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}   
};
3、算法实现
//我们已经知道，N点FFT从左到右共有log2N级蝶形，每级有N/2L组，每组有L个。
//所以FFT的C语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算（整个FFT共log2N级），
//中间层循环完成每一组的蝶形运算（每一级有N/2L组），最内层循环完成单独1个蝶形运算（每一组有L个）。
  /***【快速傅里叶变换】***/
  void FFT(void)
  {
   unsigned int i,j,k,l;
   complex top,bottom,xW;
   Reverse(); //码位倒序
   for(i=0;i<log2N;i++)   /*共log2N级*/
   {    //一级蝶形运算
    l=1<<i;//l等于2的i次方
    for(j=0;j<N;j+=2*l)  /*每L个蝶形是一组，每级有N/2L组*/
    {   //一组蝶形运算
     for(k=0;k<l;k++)   /*每组有L个*/
     {  //一个蝶形运算
       xW=mul(x[j+k+l],WN[N/(2*l)*k]); //碟间距为l
       top=add(x[j+k],xW); //每组的第k个蝶形
       bottom=sub(x[j+k],xW);
       x[j+k]=top;
       x[j+k+l]=bottom;
     }
    }
   }
  }
三、FFT计算结果验证
随便输入一个64点序列，例如
x[N]={{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0}};
在Keil中Debug查看数组变量x的FFT计算结果并与MATLAB计算结果进行比对，可以发现非常准确，说明程序编写正确！

好资料一份 下载学习一下

很好的资料，，优化得很好。仔细看了下，如果旋转因子表能用代码初始化生成就更好了移植了

参考参考，正在研究

很好很好

正在研究

Thank you! 好东西！！！！！！！！！！

好东西 收藏

为啥把文件里面的战舰V3的那个代码直接下载过去运行不了呢？求教

就是那个基于战舰V3的音乐频谱里面那个工程

感谢分享资料，有你更精彩！

真的是好资料，需要认真学习，谢谢楼主分享

感谢了 楼主

好  多料，下载学习

感谢分享资料，真的是好资料，谢谢楼主！

这个FFT的原创请私信我,我要咨询.还有懂这个的,也可以私信我,我的工作就是将你们的知识和兴趣变成钱,我们在广州,机场和高铁站边上,有工厂,免费的驿站,欢迎来坐坐.(这个FFT是要应用在STC的单片机上,做音乐频谱LED显示.公司有硬软件工程师,但他们的工作太忙,欢迎外援支持,动动键盘就行,别的有人做)

谢谢楼主分享，有学习到

学习，谢谢分享！

有没有仿真图  或者电路图

这个程序有bug,需要初始化temp=0；
第二个程序：我写的
#include <reg52.h>
unsigned char mr;
unsigned char invertir_byte (mr) {
unsigned char temp;
if(mr&0x80){temp=temp|0x01;}
if(mr&0x40){temp=temp|0x02;}
if(mr&0x20){temp=temp|0x04;}
if(mr&0x10){temp=temp|0x08;}
if(mr&0x08){temp=temp|0x10;}
if(mr&0x04){temp=temp|0x20;}
if(mr&0x02){temp=temp|0x40;}
if(mr&0x01){temp=temp|0x80;}
return (temp);
}
void main()
{
while(1)
{
P1=invertir_byte(0x33);
}
}

Thank you! 好东西！

有对比有分析，值得学习！棒棒哒！

stm32 有#include <reg52.h>？