有趣|如何用低成本MCU实现音乐频谱显示

NXP管管

有趣|如何用低成本MCU实现音乐频谱显示

♦效果展示
下面是实际频谱显示的测试效果。

如果无法播放，请点击此处在新窗口打开

最近在做一个有趣的小项目，其中有一小部分的内容的是使用FFT做音乐频谱显示。于是就有了下面这个音乐频谱显示的低成本方案，话不多说看看低成本MCU如何实现FFT音乐频谱显示吧。

♦音频采集硬件电路
音频采集的硬件电路比较简单，主要的器件就是麦克风和LM358运放。

图中电路R5可调电阻的作用是来调节运放的增益。R4的作用的是给运放一个VDD*R4/(R3+R4) 的直流偏置，这里加直流偏置是由于ADC只能采集正电压值，为了不丢失负电压的音频信号，给信号整体加了一个直流偏置。


但是这个图还有一个小问题，运放的输出端加了一个电容C2，C2会把直流偏置给隔掉。在设计时，这个电容可以去掉。


下图是按照上图搭建的音频采集电路的输出信号，图中波动信号是施加的外部音频，是我们需要做音乐频谱显示需要的信号。该信号有一个2.3v的直流偏置，在后续处理时需要减去这个偏置。

♦CTimer+ADC+DMA 音频信号采集
为了呼应标题，我们选择的MCU是LPC845，这是NXP的一款低成本的MCU。考虑到我们平常听的音乐频率大都低于5kHz，在软件设计时设置ADC采样频率为10kHz。不要问为什么，问就是采样定理。

LPC845的ADC有8个触发源，我们使用CTiimer match3来触发ADC，将寄存器SEQA_CTRL的bit 14:12设置为0x5。CTimer match 3的输出频率为10kHz。

为了确保我们采集数据的实时性，DMA建议配置成双buffer模式，以防止采样的数据被覆盖掉。
♦FFT音频信号处理
在DMA搬运ADC采样值时，使用了双buffer来搬，ADC采样值需要减去一个2.3V的直流偏置。Samples[]数组用于FFT计算。

1.    //Calculate the FFT input buffer
   
2.    if(g_DmaTransferDoneFlag_A == true)
   
3.    {
   
4.            for (i=0; i<128; i++)
   
5.            {
   
6.                   Samples[i] =(int16_t)(((g_AdcConvResult_A[i] & 0xfff0) >> 4) - 2979);//substract the 2.3v offset in the Amplifier output
   
7.             }
   
8.           g_DmaTransferDoneFlag_A = false;
   
9. 
   
10.     }
    
11.    else if(g_DmaTransferDoneFlag_B == true)
    
12.    {
    
13.            for (i=0; i<128; i++)
    
14.            {
    
15.                   Samples[i] =(int16_t)(((g_AdcConvResult_B[i] & 0xfff0) >> 4) - 2979);//substract the 2.3v offset in the Amplifier output
    
16.            }
    
17.            g_DmaTransferDoneFlag_B = false;
    
18.    }

复制代码

根据FFT算法的原理，在进行FFT计算之前，还需要将ADC的采样值Samples[]乘上一个窗函数，这里我们使用的汉宁窗函数，由于篇幅限制，具体原理可以去查看FFT算法相关的资料。

1.   //If 'Window' isn't rectangular, apply window
   
2.     if(Window == Triangular){
   
3.         //Apply a triangular window to the data.
   
4.         for(Cnt = 0; Cnt<Len; Cnt++){
   
5.             if(Cnt<(Len/2)) Samples[Cnt] = ((int32_t)Samples[Cnt]*Cnt)>>L2Len;
   
6.             else Samples[Cnt] = ((int32_t)Samples[Cnt]*((Len/2)-Cnt))>>L2Len;
   
7.         }
   
8.     }
   
9.     else if(Window == Hann){
   
10.         //Use the cosine window wavetable to apply a Hann windowing function to the samples
    
11.         for(Cnt = 0; Cnt<Len; Cnt++){
    
12.             Index = (Cnt*CWLen)>>L2Len;
    
13.             Samples[Cnt] = ((int32_t)Samples[Cnt]*(int32_t)CosWindow[Index])>>(CWBD);
    
14.         }
    
15.     }

复制代码

前面说了这么多，FFT算法才是实现音乐频谱显示的关键部分（其实上边每一步都缺一不可）。

我在网上找了好多FFT算法的资料，大家在做频谱显示时，用到最多的就是CMSIS DSP的算法库。于是乎，采用CMSIS DSP的库貌似是首选。

但是不用不知道，一用才发现，由于CMSIS DSP的库使用的是查表的方式，我的64K Flash的LPC845轻轻松松就被撑爆了。没办法，只能改用其他方案。经过不懈的查阅资料，在GitHub找到一份FFT算法的代码，这个代码写的非常简洁，而且用起来很好用，感谢发布者pyrohaz，下面是FFT代码的一部分。

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

复制代码

fix_fft(short fr[], short fi[], short m, short inverse)函数，FFT计算函数

1. int fix_fft(short fr[], short fi[], short m, short inverse)
   
2. {
   
3.     int mr, nn, i, j, l, k, istep, n, scale, shift;
   
4.     short qr, qi, tr, ti, wr, wi;
   
5. 
   
6.     n = 1 << m;
   
7. 
   
8.     /* max FFT size = N_WAVE */
   
9.     if (n > N_WAVE)
   
10.         return -1;
    
11. 
    
12.     mr = 0;
    
13.     nn = n - 1;
    
14.     scale = 0;
    
15. 
    
16.     /* decimation in time - re-order data */
    
17.     for (m=1; m<=nn; ++m) {
    
18.         l = n;
    
19.         do {
    
20.             l >>= 1;
    
21.         } while (mr+l > nn);
    
22.         mr = (mr & (l-1)) + l;
    
23. 
    
24.         if (mr <= m)
    
25.             continue;
    
26.         tr = fr[m];
    
27.         fr[m] = fr[mr];
    
28.         fr[mr] = tr;
    
29.         ti = fi[m];
    
30.         fi[m] = fi[mr];
    
31.         fi[mr] = ti;
    
32.     }      

复制代码

接 fix_fft(short fr[], short fi[], short m, short inverse)函数

1.   l = 1;
   
2.     k = LOG2_N_WAVE-1;
   
3.     while (l < n) {
   
4.         if (inverse) {
   
5.             /* variable scaling, depending upon data */
   
6.             shift = 0;
   
7.             for (i=0; i<n; ++i) {
   
8.                 j = fr[i];
   
9.                 if (j < 0)
   
10.                     j = -j;
    
11.                 m = fi[i];
    
12.                 if (m < 0)
    
13.                     m = -m;
    
14.                 if (j > 16383 || m > 16383) {
    
15.                     shift = 1;
    
16.                     break;
    
17.                 }
    
18.             }
    
19.             if (shift)
    
20.                 ++scale;
    
21.         } else {
    
22.             /*
    
23.                 fixed scaling, for proper normalization --
    
24.                 there will be log2(n) passes, so this results
    
25.                 in an overall factor of 1/n, distributed to
    
26.                 maximize arithmetic accuracy.
    
27.             */
    
28.             shift = 1;
    
29.         }

复制代码

接fix_fftr(short f[], int m, int inverse)函数

1.   /*
   
2.             it may not be obvious, but the shift will be
   
3.             performed on each data point exactly once,
   
4.             during this pass.
   
5.         */
   
6.         istep = l << 1;
   
7.         for (m=0; m<l; ++m) {
   
8.             j = m << k;
   
9.             /* 0 <= j < N_WAVE/2 */
   
10.             wr =  Sinewave[j+N_WAVE/4];
    
11.             wi = -Sinewave[j];
    
12.             if (inverse)
    
13.                 wi = -wi;
    
14.             if (shift) {
    
15.                 wr >>= 1;
    
16.                 wi >>= 1;
    
17.             }
    
18.             for (i=m; i<n; i+=istep) {
    
19.                 j = i + l;
    
20.                 tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
    
21.                 ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
    
22.                 qr = fr[i];
    
23.                 qi = fi[i];
    
24.                 if (shift) {
    
25.                     qr >>= 1;
    
26.                     qi >>= 1;
    
27.                 }
    
28.                 fr[j] = qr - tr;
    
29.                 fi[j] = qi - ti;
    
30.                 fr[i] = qr + tr;
    
31.                 fi[i] = qi + ti;
    
32.             }
    
33.         }
    
34.         --k;
    
35.         l = istep;
    
36.     }
    
37.     return scale;
    
38. }

复制代码

1. /*
   
2.   fix_fftr() - forward/inverse FFT on array of real numbers.
   
3.   Real FFT/iFFT using half-size complex FFT by distributing
   
4.   even/odd samples into real/imaginary arrays respectively.
   
5.   In order to save data space (i.e. to avoid two arrays, one
   
6.   for real, one for imaginary samples), we proceed in the
   
7.   following two steps: a) samples are rearranged in the real
   
8.   array so that all even samples are in places 0-(N/2-1) and
   
9.   all imaginary samples in places (N/2)-(N-1), and b) fix_fft
   
10.   is called with fr and fi pointing to index 0 and index N/2
    
11.   respectively in the original array. The above guarantees
    
12.   that fix_fft "sees" consecutive real samples as alternating
    
13.   real and imaginary samples in the complex array.
    
14. */
    
15. int fix_fftr(short f[], int m, int inverse)
    
16. {
    
17.     int i, N = 1<<(m-1), scale = 0;
    
18.     short tt, *fr=f, *fi=&f[N];
    
19. 
    
20.     if (inverse)
    
21.         scale = fix_fft(fi, fr, m-1, inverse);
    
22.     for (i=1; i<N; i+=2) {
    
23.         tt = f[N+i-1];
    
24.         f[N+i-1] = f[i];
    
25.         f[i] = tt;
    
26.     }
    
27.     if (! inverse)
    
28.         scale = fix_fft(fi, fr, m-1, inverse);
    
29.     return scale;
    
30. }

复制代码

int fix_fft(short fr[], short fi[], short m, short inverse) 是FFT算法的计算函数，fr[]是ADC采集到信号值的实部，fi[]是ADC采集到信号值的虚部。经过fix_fft函数处理之后，fr[]是FFT计算所得实部，fi[]是计算所得的虚部。

我们最终要显示的音乐频谱其实是FFT频域中音频的幅值，幅值的计算是实部的平方+虚部的平方开根号。下面是具体的幅值计算部分的代码，每一个幅值点对应OLED的一列像素点。

1. //Calculate the magnitude
   
2.     for(Cnt = 0; Cnt<Len; Cnt++){
   
3.         //Square the real components
   
4.         R = (int32_t)Samples[Cnt]*(int32_t)Samples[Cnt];
   
5. 
   
6.         //Square the imaginary components
   
7.         I = (int32_t)Imag[Cnt]*(int32_t)Imag[Cnt];
   
8. 
   
9.         //Calculate the magnitude and sum into BufSum. As there
   
10.         //are more samples than columns, sum the magnitude of
    
11.         //the bins. E.g. for 128 samples and 64 pixels on the LCD,
    
12.         //each column should be the sum of 2 bins (128/64). If there
    
13.         //are 512 samples and 64 pixels on the LCD, each column should
    
14.         //be the sum of 8 bins (512/64)
    
15.         BufSum += isqrt(R + I);
    
16. 
    
17.         //Calculate the index for the column
    
18.         Index = Cnt;
    
19. 
    
20.         //Reset BufSum once the index has changed
    
21.         if(Index != IndO){
    
22.             //Write buffer sum to the column. Scale BufSum dependent on
    
23.             //input amplitude range. If 'ColumnFilter' is more than 0,
    
24.             //low pass filter the value to smooth out erratic changes.
    
25.             if(Decibels){
    
26.                 Col[Index] += (int32_t)(20.0f*log10f(BufSum)-Col[Index])>>ColumnFilter; //calculate the DB
    
27.             }
    
28.             else{
    
29.                 Col[Index] += (BufSum-Col[Index])>>ColumnFilter;  //calculate the amplitude
    
30.             }
    
31. 
    
32.             //Limit maximum column value
    
33.             if(Col[Index] >= YPix-9) Col[Index] = YPix-10;
    
34. 
    
35.             IndO = Index;
    
36.             BufSum = 0;
    
37.         }
    
38.     }

复制代码

阅读全文

y369369

这个看着不错

流水源

有机会试试

stm1024

你们能看清到视频嘛？

NXP管管

stm1024 发表于 2022-9-22 10:57
你们能看清到视频嘛？

可以吧

wudianjun2001

很不错。。。

jundao721

NXP管管 发表于 2022-9-22 12:39
可以吧

换个计算机试试，看不到视频，漆黑一片。

sumoon_yao

视频看不到

超級稻草人

收藏了。。

lospring

做一个成品当转盘一等奖