Como Utilizar a FFT com o PIC32MX - DSP Library

Neste post você encontra informações de como utilizar a Fast Fourier Transformer (FFT) em seu projeto com PIC32MX.

A família de microcontroladores PIC32MX não possui instruções de DSP. Operações de DSP estão disponíveis na família PIC32MZ. Aqui iremos apresentar a utilização da função mips_fft16 ( ) disponível na biblioteca DSP Library para o PIC32MX. Utilizaremos a linguagem C com o compilador MPLAB XC32 no ambiente de desenvolvimento MPLAB X.

Não entraremos em detalhes da teoria da Transformada Rápida de Fourier - FFT. Partiremos direto para a utilização da função mips_fft16 ( ). Para saber mais sobre a FFT você pode acessar estes posts: "Como Utilizar a FFT com o dsPIC - DSP Library" e "FFT-PIC - Um Analisador de Espectro para a Rede de Energia Elétrica ".

No início do nosso projeto devemos incluir os headers "dsplib_dsp.h" e "fftc.h" da biblioteca DSP, conforme mostra o Código 1. O arquivo "fftc.h" contém os coeficientes utilizados nos cálculos da FFT.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <p32xxxx.h>
#include <xc.h>         // comentar esta linha para usar o MPLAB C32
#include <plib.h>

#include "dsplib_dsp.h"
#include "fftc.h"

Código 1 - Incluindo headers da DSP Library ao projeto.

Para o nosso projeto temos que definir o número de "pontos" da FFT, conforme o Código 2.

/*         FFT com 128 pontos       */

#define N    128        // N = 2^7 = 128

#define fftc fft16c128  // from fftc.h, for N = 128

Código 2 - Define 128 pontos da FFT.

128 amostras no domínio do tempo (lidas pelo conversor analógico-digital, por exemplo), se transformarão em 64 pontos no domínio da frequência.

Dados da FFT deste exemplo:
- Frequência de amostragem (Fs): 7,68 ksps;
- Pontos da FFT: 128;
- log2 (N) = 7;
- Passo entre as frequências (Fo): 60 Hz.

A frequência de amostragem (Fs) do sinal de entrada irá determinar o espaçamento entre as frequências de cada ponto de saída (Fo).

A frequência de amostragem (Fs) do nosso sinal gerado (a onda quadrada) é de 7680 Hz, portanto Fo é 60 Hz, ou seja, 7680 sps / 128 pontos da FFT!

Isto quer dizer que temos na saída o valor do espectro das frequências múltiplas de 60 Hz, de 0 à 3840 Hz (60 Hz * 64 pontos no domínio da frequência).

Para o sinal de entrada que será aplicado a FFT iremos utilizar um vetor com dados de uma onda quadrada de 780 Hz e frequência de amostragem 7680 Hz. Esta onda foi gerada com o auxílio do Software dsPICworks, da Microchip. Para mais informações de como gerar este vetor de dados acesse "Como Utilizar a FFT com o dsPIC - DSP Library". O vetor sinal_Entrada [ ] do tipo int16, com 128 amostragens é visto no Código 3.

int16 sinal_Entrada [128]= {                   //vetor que ira armazenar os dados de entrada
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,        //valores iniciais declarados para simulacao
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,        //e teste do programa
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,        //square signal 780 Hz, 7680 sps, -1.0 a +1.0
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
0x8001, 0x8001, 0x8001};

 

Código 3 - Vetor com dados de entrada para a simulação e cálculo da FFT.

 A seguir precisamos declarar os vetores e variáveis utilizadas para o cálculo da função mips_fft16 ( ), conforme o Código 4.

    int log2N = 7;      // log2(128) = 7        N = 128 = (1<<log2N)     N = 2^7 = 128
    int16c din[N];
    int16c dout[N];
    int16c scratch[N];
    int16c *p_fftc = &fftc;

Código 4 - Declaração das variáveis utilizadas no cálculo da FFT.

Os vetores utilizados pela função mips_fft16 ( ) são do tipo int16c, ou seja, inteiro 16-bit complexo, com partes real e imaginária.

No início da função main ( ) iremos transferir os coeficientes para o cálculo da FFT, conhecidos como Twiddle Factors, da memória de programa para a memória RAM, para execução mais rápida do cálculo. Veja o Código 5.

    p_fftc = &fftc;                      //inicializa ponteiro

    for (i = 0; i < N / 2; i++)          //move twiddleFactors da flash para a RAM
    {
        twiddleFactors [i] = *p_fftc;
        *p_fftc++;                       //incrementa o ponteiro
    }

Código 5 - Movendo os coeficientes para a memória RAM.

Agora podemos partir para um loop while onde continuamente o PIC32 irá mover os dados contidos no vetor sinal_Entrada [ ] para o vetor din [ ] e limitar a escala entre -0.5 e +0.5. Note que sinal_Entrada [ ] pode conter dados vindos de um conversor analógico-digital ADC, por exemplo, e din [ ] é o vetor com dados complexos utilizados nos cálculos da FFT.

   while (TRUE) {

        //move buffer do sinal de entrada (real)para vetor da FFT (complexo)
        for (i = 0; i < N; i++)
        {
            din[i].re = sinal_Entrada [i];  //parte real
            din[i].im = 0;                  //parte imaginaria
            din[i].re = din[i].re >> 1;     //divide pela metade, escalona para -0.5 a +0.5
        }
          
        ...

  

Código 6 - Preparando o vetor de dados complexos de entrada para os cálculos da FFT.

A chamada para a função mips_fft16 ( ), que aplica a Transformada Rápida de Fourier ao sinal, é mostrada no Código 7. Detalhes sobre os argumentos desta função você encontra em "PIC32 DSP Library".

        // Executa o Calculo da FFT
        mips_fft16 (dout, din, twiddleFactors, scratch, log2N);

Código 7 - Cálculo da FFT.


O resultado da FFT está contido no vetor complexo dout [ ].

Para exibir o resultado gráfico da FFT comumente se computa o single sided FFT, a magnitude quadrática do espectro de frequência, ou o módulo do número complexo. O Código 8 apresenta as fórmulas que devem ser aplicadas a cada posição do vetor de dados.

        /*     Computa o single sided FFT    */
        // Calcula a soma do quadrado das partes real e imaginária transformando vetor complexo em um vetor real.
  
        square_Magnitude = parte_Real^2 + parte_Imaginaria^2
        

        /*          Calcula o Modulo         */
        // Extrai a raiz quadrada da soma dos quadrados da parte real e imaginaria
        
        modulo = sqrtf(parte_Real^2 + parte_Imaginaria^2)  

        
  

Código 8 - Cálculo da magnitude quadrática do resultado da FFT.

Por fim, podemos calcular o espectro da frequência de maior magnitude do sinal amostrado, conforme exemplificado no Código 9.

  
       // Encontra a frequency Bin ( indice da frequencia) que possui a maior energia
        picoFrequenciaBin = vectorMaxBin (PONTOS_FFT/2,squareMagnitude);

        // Calcula a frequencia em Hz com maior espectro
        picoFrequencia = picoFrequenciaBin*(SAMPLING_FREQ/PONTOS_FFT);

Código 9 - Cálculo do espectro da frequência de maior magnitude.

O programa agora pode ser testado.


Para confirmar se o seu projeto está correto, verifique se a variável picoFrequencia apresenta o valor de 780 [Hz] após os cálculos. Esta é a frequência da onda quadrada gerada no dsPICwork para a simulação dos cálculos da FFT.

A posição do vetor com o espectro de frequência de maior valor é SquareMagnitude [13].

Veja os cálculos finais:

Fs               /    PONTOS_FFT              =      Fo
7680 sps    /     128 pontos da FFT         =     60 Hz

60 Hz * [13]       =        780 Hz!


O Vídeo 1 mostra este código rodando no MPLAB X em modo simulação passo-a-passo.

Vídeo 1 - Simulação dos cálculos da FFT com o PIC32MX.

O Vídeo 2 mostra uma aplicação da FFT com o PIC32MX para a identificação da frequência de um sinal. O display gráfico mostra os dados do vetor squareMagnitude [ ] como resultado da FFT, o espectro das frequências.

Vídeo 2 - Aplicação da FFT para identificar a frequência de um sinal.

Para quem gosta de música e áudio, o Vídeo 3 mostra uma aplicação bem legal para a FFT: um Analisador de Espectro de Áudio para a animação de um MP3 Player com o PIC32MX!

 Vídeo 3 - PIC32 MP3 Player & Audio Spectrum Analyser com a FFT.

Foi utilizada a função mips_fft16 para o cálculo da FFT com 128 pontos.

A FFT é testada com o arquivo de áudio Sine3kHz.MP3 contendo uma onda senoidal variando de 20 Hz à 3 kHz.

Referências:

HAYKIN, Simon S.; VAN VEEN, Barry. Sinais e sistemas. Porto Alegre: Bookman, 2001. 668p.

MICROCHIP. Code Examples. CE018 - Using the Fast Fourier Transform (FFT) for Frequency Detection. Disponível em: http://ww1.microchip.com/downloads/en/DeviceDoc/CE018_FFT_DSPlib_101007.zip

MICROCHIP. Software Libraries. Disponível em: http://www.microchip.com/SoftwareLib.aspx

MICROCHIP. PIC32 DSP Library. Disponível em: http://www.microchip.com/stellent/idcplg?IdcService=SS_GET_PAGE&nodeId=2680&dDocName=en552827