音響工程—FFT頻譜

ID:75263 · 發表于 2015-6-9 04:01

大學的第一次電子設計競賽初賽的題目就是《音頻頻譜柱狀顯示電路》，當時貪玩，沒能按時完成作品，因此連晉級決賽的名額都沒有�，F在回想起來，真的覺得自己很不爭氣。

圖1 2011年電子設計競賽初賽題目

第二年老師重新布置任務，我接到的還是繼續完成這個題目，硬件部分的濾波器倒是做好了，并且龍富做好了后期的顯示條部分，但因期末考試的原因沒能把硬件結合起來。
前些天翻了一下以前的筆記本，一幕幕記憶彌漫出來。剛好有點時間，整理一下思路，決定用FFT算法去實現這個功能。 FFT是離散傅立葉變換的快速算法，可以將一個信號變換到頻域。有些信號在時域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。想想《高數》、《復變函數》、《信號系統》、《現代自動控制系統》、《離散自動控制》都有涉及到，也算是專業知識了，今天就用它來實現多通道帶通濾波器的數據解析。
本設計是基于單片機C的FFT快速傅里葉變換，用于音頻信號分析的一個典型音頻頻譜顯示儀。整體實現框圖如下：

圖2 整體框架

MCU控制部分使用32位ARM單片機，7寸TFT彩屏顯示器，FFT算法使用純C代碼實現。整個實現過程中關鍵的一部分是數據的采樣，為了完整的采樣到32768HZ的音頻信號，采用定時器觸發ADC，并使用DMA傳送到指定內存區域，再進行FFT運算。幾天的FFT學習和數據分析，終于理解了整個功能的實現過程。
關鍵部分的是FFT核心轉換程序，整理如下（聲明：代碼摘自網上）：
FFT.h頭文件
//                      快速福利葉變換C函數
//函數簡介：此函數是通用的快速傅里葉變換C語言函數，移植性強，以下部分不依
//       賴硬件。此函數采用聯合體的形式表示一個復數，輸入為自然順序的復
//       數（輸入實數是可令復數虛部為0），輸出為經過FFT變換的自然順序的
//       復數
//使用說明：使用此函數只需更改宏定義FFT_N的值即可實現點數的改變，FFT_N的
//       應該為2的N次方，不滿足此條件時應在后面補0
//函數調用：FFT(s);
//時間：2010-2-20
//版本：Ver1.0
//參考文獻：
#ifndef __FFT_H
#define __FFT_H
#include "math.h"
#define PI 3.1415926535897932384626433832795028841971//定義圓周率值
#define FFT_N512                                  //定義福利葉變換的點數
struct compx {floatreal,imag;};                //定義一個復數結構
extern struct compxs[FFT_N];                   //FFT輸入和輸出：從S[1]開始存放，根據大小自己定義
void FFT(struct compx*xin);                   //FFT核心算法
#endif

FFT.c文件：
#include "FFT.h"
struct compx s[FFT_N];//FFT輸入和輸出：從S[1]開始存放，根據大小自己定義

struct compx EE(struct compx a,struct compxb)
{
struct compx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}

void FFT(struct compx *xin)
{
  int f,m,nv2,nm1,i,k,l,j=0;
  struct compx u,w,t;

  nv2=FFT_N/2;                //變址運算，即把自然順序變成倒位序，采用雷德算法
  nm1=FFT_N-1;
for(i=0;i
{
if(i
{
   t=xin[j];
   xin[j]=xin[ i];
   xin[ i]=t;
}
k=nv2;                         //求j的下一個倒位序
while(k<=j)                //如果k<=j,表示j的最高位為1
{
   j=j-k;                         //把最高位變成0
   k=k/2;                      //k/2，比較次高位，依次類推，逐個比較，直到某個位為0
}
  j=j+k;                         //把0改為1
  }

  {
intle,lei,ip;                   //FFT運算核，使用蝶形運算完成FFT運算
f=FFT_N;
  for(l=1;(f=f/2)!=1;l++)  //計算l的值，即計算蝶形級數
      ;
for(m=1;m<=l;m++)    // 控制蝶形結級數
  {                                     //m表示第m級蝶形，l為蝶形級總數l=log（2）N
le=2<<(m-1);             //le蝶形結距離，即第m級蝶形的蝶形結相距le點
lei=le/2;                      //同一蝶形結中參加運算的兩點的距離
u.real=1.0;                   //u為蝶形結運算系數，初始值為1
u.imag=0.0;
w.real=cos(PI/lei);       //w為系數商，即當前系數與前一個系數的商
w.imag=-sin(PI/lei);
for(j=0;j<=lei-1;j++)    //控制計算不同種蝶形結，即計算系數不同的蝶形結
{
   for(i=j;i<=FFT_N-1;i=i+le)  //控制同一蝶形結運算，即計算系數相同蝶形結
   {
   ip=i+lei;                   //i，ip分別表示參加蝶形運算的兩個節點
   t=EE(xin[ip],u);          //蝶形運算，詳見公式
   xin[ip].real=xin[ i].real-t.real;
   xin[ip].imag=xin[ i].imag-t.imag;
   xin[ i].real=xin[ i].real+t.real;
   xin[ i].imag=xin[ i].imag+t.imag;
   }
   u=EE(u,w);                   //改變系數，進行下一個蝶形運算
}
}
  }
}

FFT分析后的數據雖然有一定的規律，但是還不能滿足頻譜顯示的需求，需在程序控制上針對顯示器的尺寸進行數據放縮，頂值限幅等才能輸出酷炫的FFT頻譜界面。
汗水終究交換了成功，FFT頻譜顯示，如圖：