MCU如何實(shí)現(xiàn)FFT音樂頻譜顯示

最近在做一個有趣的小項(xiàng)目，其中有一小部分的內(nèi)容的是使用FFT做音樂頻譜顯示。于是就有了下面這個音樂頻譜顯示的低成本方案，話不多說看看低成本MCU如何實(shí)現(xiàn)FFT音樂頻譜顯示吧。

音頻采集硬件電路

音頻采集的硬件電路比較簡單，主要的器件就是麥克風(fēng)和LM358運(yùn)放。

圖中電路R5可調(diào)電阻的作用是來調(diào)節(jié)運(yùn)放的增益。R4的作用的是給運(yùn)放一個VDD*R4/(R3+R4) 的直流偏置，這里加直流偏置是由于ADC只能采集正電壓值，為了不丟失負(fù)電壓的音頻信號，給信號整體加了一個直流偏置。

但是這個圖還有一個小問題，運(yùn)放的輸出端加了一個電容C2，C2會把直流偏置給隔掉。在設(shè)計時，這個電容可以去掉。

下圖是按照上圖搭建的音頻采集電路的輸出信號，圖中波動信號是施加的外部音頻，是我們需要做音樂頻譜顯示需要的信號。該信號有一個2.3v的直流偏置，在后續(xù)處理時需要減去這個偏置。

CTimer+ADC+DMA 音頻信號采集

為了呼應(yīng)標(biāo)題，我們選擇的MCU是LPC845，這是NXP的一款低成本的MCU?？紤]到我們平常聽的音樂頻率大都低于5kHz，在軟件設(shè)計時設(shè)置ADC采樣頻率為10kHz。不要問為什么，問就是采樣定理。

LPC845的ADC有8個觸發(fā)源，我們使用CTiimer match3來觸發(fā)ADC，將寄存器SEQA_CTRL的bit 14:12設(shè)置為0x5。CTimer match 3的輸出頻率為10kHz。

為了確保我們采集數(shù)據(jù)的實(shí)時性，DMA建議配置成雙buffer模式，以防止采樣的數(shù)據(jù)被覆蓋掉。

FFT音頻信號處理

在DMA搬運(yùn)ADC采樣值時，使用了雙buffer來搬，ADC采樣值需要減去一個2.3V的直流偏置。Samples[]數(shù)組用于FFT計算。

    //Calculate the FFT input buffer
   if(g_DmaTransferDoneFlag_A == true)
   {
   for (i=0; i<128; i++)
   {
  Samples[i] =(int16_t)(((g_AdcConvResult_A[i] & 0xfff0) >> 4) - 2979);//substract the 2.3v offset in the Amplifier output
    }
  g_DmaTransferDoneFlag_A = false;

    }
   else if(g_DmaTransferDoneFlag_B == true)
   {
   for (i=0; i<128; i++)
   {
  Samples[i] =(int16_t)(((g_AdcConvResult_B[i] & 0xfff0) >> 4) - 2979);//substract the 2.3v offset in the Amplifier output
   }
   g_DmaTransferDoneFlag_B = false;
   }

根據(jù)FFT算法的原理，在進(jìn)行FFT計算之前，還需要將ADC的采樣值Samples[]乘上一個窗函數(shù)，這里我們使用的漢寧窗函數(shù)，由于篇幅限制，具體原理可以去查看FFT算法相關(guān)的資料。

    //If 'Window' isn't rectangular, apply window
    if(Window == Triangular){
        //Apply a triangular window to the data.
        for(Cnt = 0; Cnt>L2Len;
            else Samples[Cnt] = ((int32_t)Samples[Cnt]*((Len/2)-Cnt))>>L2Len;
        }
    }
    else if(Window == Hann){
        //Use the cosine window wavetable to apply a Hann windowing function to the samples
        for(Cnt = 0; Cnt>L2Len;
            Samples[Cnt] = ((int32_t)Samples[Cnt]*(int32_t)CosWindow[Index])>>(CWBD);
        }
    }

前面說了這么多，F(xiàn)FT算法才是實(shí)現(xiàn)音樂頻譜顯示的關(guān)鍵部分（其實(shí)上邊每一步都缺一不可）。

我在網(wǎng)上找了好多FFT算法的資料，大家在做頻譜顯示時，用到最多的就是CMSIS DSP的算法庫。于是乎，采用CMSIS DSP的庫貌似是首選。

但是不用不知道，一用才發(fā)現(xiàn)，由于CMSIS DSP的庫使用的是查表的方式，我的64K Flash的LPC845輕輕松松就被撐爆了。沒辦法，只能改用其他方案。經(jīng)過不懈的查閱資料，在GitHub找到一份FFT算法的代碼，這個代碼寫的非常簡潔，而且用起來很好用，感謝發(fā)布者pyrohaz，下面是FFT代碼的一部分。

/*
  FIX_MPY() - fixed-point multiplication & scaling.
  Substitute inline assembly for hardware-specific
  optimization suited to a particluar DSP processor.
  Scaling ensures that result remains 16-bit.
*/
inline short FIX_MPY(short a, short b)
{
    /* shift right one less bit (i.e. 15-1) */
    int c = ((int)a * (int)b) >> 14;
    /* last bit shifted out = rounding-bit */
    b = c & 0x01;
    /* last shift + rounding bit */
    a = (c >> 1) + b;
    return a;
}

fix_fft(short fr[], short fi[], short m, short inverse)函數(shù)，F(xiàn)FT計算函數(shù)

int fix_fft(short fr[], short fi[], short m, short inverse)
{
    int mr, nn, i, j, l, k, istep, n, scale, shift;
    short qr, qi, tr, ti, wr, wi;

    n = 1 << m;

    /* max FFT size = N_WAVE */
    if (n > N_WAVE)
        return -1;

    mr = 0;
    nn = n - 1;
    scale = 0;

    /* decimation in time - re-order data */
    for (m=1; m<=nn; ++m) {
        l = n;
        do {
            l >>= 1;
        } while (mr+l > nn);
        mr = (mr & (l-1)) + l;

        if (mr <= m)
            continue;
        tr = fr[m];
        fr[m] = fr[mr];
        fr[mr] = tr;
        ti = fi[m];
        fi[m] = fi[mr];
        fi[mr] = ti;
    }       

接 fix_fft(short fr[], short fi[], short m, short inverse)函數(shù)

    l = 1;
    k = LOG2_N_WAVE-1;
    while (l < n) {
        if (inverse) {
            /* variable scaling, depending upon data */
            shift = 0;
            for (i=0; i 16383 || m > 16383) {
                    shift = 1;
                    break;
                }
            }
            if (shift)
                ++scale;
        } else {
            /*
                fixed scaling, for proper normalization --
                there will be log2(n) passes, so this results
                in an overall factor of 1/n, distributed to
                maximize arithmetic accuracy.
            */
            shift = 1;
        }

接fix_fftr(short f[], int m, int inverse)函數(shù)

        /*
            it may not be obvious, but the shift will be
            performed on each data point exactly once,
            during this pass.
        */
        istep = l << 1;
        for (m=0; m>= 1;
                wi >>= 1;
            }
            for (i=m; i>= 1;
                    qi >>= 1;
                }
                fr[j] = qr - tr;
                fi[j] = qi - ti;
                fr[i] = qr + tr;
                fi[i] = qi + ti;
            }
        }
        --k;
        l = istep;
    }
    return scale;
}

/*
  fix_fftr() - forward/inverse FFT on array of real numbers.
  Real FFT/iFFT using half-size complex FFT by distributing
  even/odd samples into real/imaginary arrays respectively.
  In order to save data space (i.e. to avoid two arrays, one
  for real, one for imaginary samples), we proceed in the
  following two steps: a) samples are rearranged in the real
  array so that all even samples are in places 0-(N/2-1) and
  all imaginary samples in places (N/2)-(N-1), and b) fix_fft
  is called with fr and fi pointing to index 0 and index N/2
  respectively in the original array. The above guarantees
  that fix_fft "sees" consecutive real samples as alternating
  real and imaginary samples in the complex array.
*/
int fix_fftr(short f[], int m, int inverse)
{
    int i, N = 1<<(m-1), scale = 0;
    short tt, *fr=f, *fi=&f[N];

    if (inverse)
        scale = fix_fft(fi, fr, m-1, inverse);
    for (i=1; i

	

	int fix_fft(short fr[], short fi[], short m, short inverse) 是FFT算法的計算函數(shù)，fr[]是ADC采集到信號值的實(shí)部，fi[]是ADC采集到信號值的虛部。經(jīng)過fix_fft函數(shù)處理之后，fr[]是FFT計算所得實(shí)部，fi[]是計算所得的虛部。

	我們最終要顯示的音樂頻譜其實(shí)是FFT頻域中音頻的幅值，幅值的計算是實(shí)部的平方+虛部的平方開根號。下面是具體的幅值計算部分的代碼，每一個幅值點(diǎn)對應(yīng)OLED的一列像素點(diǎn)。

	
    //Calculate the magnitude
    for(Cnt = 0; Cnt>ColumnFilter; //calculate the DB
            }
            else{
                Col[Index] += (BufSum-Col[Index])>>ColumnFilter;  //calculate the amplitude
            }

            //Limit maximum column value
            if(Col[Index] >= YPix-9) Col[Index] = YPix-10;

            IndO = Index;
            BufSum = 0;
        }
    }


	

	

	   審核編輯：彭靜