Mostrando postagens com marcador FFT. Mostrar todas as postagens
Mostrando postagens com marcador FFT. Mostrar todas as postagens

terça-feira, 12 de maio de 2015

Como utilizar IIR Filter com o PIC32MX - DSP Library

Este artigo demonstra como elaborar um Filtro Digital IIR utilizando o microcontrolador PIC32MX com a DSP Library da Microchip.

Filtros Digitais do tipo FIR (Resposta ao Impulso Finita) e IIR (Resposta ao Impulso Infinita) são muito utilizados por DSPs (Processadores Digitais de Sinais) para a filtragem de sinais.

Semelhantes aos filtros analógicos, os filtros digitais podem ser classificados em 4 tipos quanto à frequência:
- Passa-baixa (lowpass);
- Passa-alta (highpass);
- Passa-banda (bandpass);
- Rejeita-banda (bandstop).

Funções de Janela são utilizadas em Filtros Digitais FIR podendo ser:
- janela Retangular;
- janela de Hamming;
- janela de Blackman;
- janela de Kaiser.

Filtros Digitais IIR podem utilizar modelagem dos filtros analógicos:
- filtro Butterworth
- filtro Tschebyscheff
- filtro Tschebyscheff inverso
- filtro Elliptic
- filtro Bessel


De um modo geral os Filtros Digitais IIR são representados pela Equação 1.


Equação 1 - Equação básica para um Filtro Digital.


 A Função de Transferência H (z) para um Filtro Digital é vista na Equação 2.

 Equação 2 - Função de Transferência H (z) de um Filtro Digital. Fonte: ORFANIDS.



Considerando a0 =1  e M=N=2 obtemos a Equação 3 para um FIltro Digital IIR de Segunda Ordem.

Equação 3 - Função de Transferência para um Filtro Digital de Segunda Ordem. Fonte: ORFANIDS.


A Equação 4 mostra a equação diferencial para o FIltro Digital IIR de Segunda Ordem ou Biquadrado.

Equação 4 - Equação Diferencial para o Filtro Digital de Segunda Ordem. Fonte: ORFANIDS.


A equação diferencial do IIR FIlter Biquad é representada pelo Fluxograma 1.

Fluxograma 1 - Realização do Filtro IIR Segunda Ordem na Forma Direta. Fonte: ORFANIDS


Cada seção biquadrada do FIltro IIR da DSP Library da Microchip para o PIC32MX opera segundo o Fluxograma 2 (Second  Order Direct Form II Structure Transposed).
Fluxograma 2 - Operação do FIltro IIR da DSP Library da Microchip. Fonte: Microchip.



Criando um Filtro Digital IIR


A Figura 1 mostra os 3 passos para você criar seu Filtro Digital.



 Figura 1 - Passos para criar um Filtro Digital.

Este artigo mostra um exemplo para a construção de 4 tipos de Filtros Digitais IIR:
a) filtro passa-baixa de 1000 Hz;
b) filtro passa-alta de 2000 Hz;
c) filtro passa-banda de 1500 Hz à 2000 Hz;
d) filtro rejeita-banda de 2000 Hz à 2500 Hz.


Vamos a algumas considerações do nosso projeto.

Objetivo do circuito: 

Filtrar um sinal entre 0 e 3840 Hz. Aplicar os 4 tipos de filtros IIR para verificar a resposta de cada um deles.

Especificações: 

- o sinal de entrada é uma onda senoidal gerada pelo próprio microcontrolador através da equação trigonométrica  cos(a+b) = 2 cos a cos b - cos(a-b).

- de acordo com o Teorema de Nyquist, a frequência de amostragem do sinal de entrada será 7680 Hz.

- o microcontrolador escolhido é o PIC32MX, um microcontrolador de 32-bit da Microchip, que possui as funções para o Filtro Digital IIR disponíveis na DSP Library. Possui 80 MHz, 512 kB de memória Flash e 128 KB High Speed SRAM.

- será utilizado o dsPIC Filter Design Lite, um programa de Digital Filter Design, para definir os parâmetros do filtro digital IIR a fim de obter o sinal filtrado especificado. O software fornece os coeficientes necessários para o cálculo do filtro.

- os coeficientes são utilizados no código fonte C do PIC32. O programa consiste basicamente em declarar as variáveis, inicializar os ponteiros e variáveis, iniciar a estrutura do filtro com mips_iir16_setup  ( ) e executar o cálculo do filtro digital com mips_iir16 ( ).

- é aplicada a transformada Rápida de Fourier (FFT) ao sinal filtrado para verificar a resposta do filtro.

- para teste do nosso circuito o PIC32 varia a frequência do sinal gerado de 0 à 3840 Hz.

- o sinal senoidal gerado e filtrado assim como a FFT são mostrados em um display gráfico.


Digital Filter Design: definindo os parâmetros do Filtro Digital

No primeiro passo temos que definir alguns parâmetros para poder utilizar um software capaz de gerar os coeficientes necessários para o cálculo do filtro digital pelo PIC.

Utilizaremos o Filtro Digital IIR Biquadrado (cascata de segunda ordem). Cada IIR FIlter Biquad possui 2 pólos e 2 zeros.

O FIltro Analógico utilizado foi o Butterworth de 4ª ordem. 2 filtros de segunda ordem (biquadrado) em cascata. No total teremos 4 pólos e 4 zeros.
 
Veja a seguir o Design elaborado para os 4 tipos de Filtros utilizando o dsPIC Filter Design Lite: lowpass, highpass, bandpass e bandstop.

 Lowpass IIR FIlter Design


 A Figura 2 traz graficamente a representação de cada parâmetro do filtro:
- frequência de amostragem;
- passband;
- stopband;
- máximo ripple no passband;
- mínima atenuação do stopband. 


                                 Figura 2 - Especificação do FIltro Digital passa-baixa.

No segundo passo utilizaremos o dsPIC FD Lite, o Digital Filter Design Software,  para modelar o nosso filtro digital com estes parâmetros e gerar os coeficientes para ser utilizado no programa do PIC32.

O método de realização do filtro, quantização, neste caso é "cascated second order sections". Estamos projetando um IIR FIlter Biquad, isto é, com 2 pólos e 2 zeros para cada seção do fiiltro.

No menu Design / IIR Design... iniciaremos o design para um Filtro Digital IIR. 

Escolheremos o tipo de filtro: lowpass.

Na segunda tela, Figura 3, preencheremos os parâmetros do filtro passa-baixa conforme ilustrado anteriormente na Figura 2.




Figura 3 - Configurando os parâmetros para o filtro passa-baixa.

A Figura 4 mostra o resultado matemático para a modelagem do nosso filtro digital IIR passa-baixa de acordo com os parâmetros fornecidos.

Figura 4 - Resultado da Modelagem do Filtro Digital IIR lowpass.


 Highpass IIR FIlter Design

A Figura 5 traz graficamente a representação de cada parâmetro do filtro highpass.


Figura 5 - Especificação do FIltro Digital passa-alta.

Os parâmetros do filtro passa-alta são configurados conforme a Figura 6.
 
Figura 6 - Configurando os parâmetros para o filtro passa-alta.

 A Figura 7 mostra o resultado matemático para a modelagem do filtro digital IIR passa-alta.

Figura 7 - Resultado da Modelagem do Filtro Digital IIR highpass.


 Bandpass IIR FIlter Design

 A Figura 8 traz graficamente a representação de cada parâmetro do filtro bandpass.


Figura 8 - Especificação do FIltro Digital passa-banda.

Os parâmetros do filtro passa-banda são configurados conforme a Figura 9.
 
Figura 9 - Configurando os parâmetros para o filtro passa-banda.

 A Figura 10 mostra o resultado matemático para a modelagem do filtro digital IIR passa-banda.

Figura 10 - Resultado da Modelagem do Filtro Digital IIR bandpass.


 Bandstop IIR FIlter Design

 A Figura 11 traz graficamente a representação de cada parâmetro do filtro bandstop.

Figura 11 - Especificação do FIltro Digital rejeita-banda.

Os parâmetros do filtro rejeita-banda são configurados conforme a Figura 12.


Figura 12 - Configurando os parâmetros para o filtro rejeita-banda.

 A Figura 13 mostra o resultado matemático para a modelagem do filtro digital IIR rejeita-banda.

Figura 13 - Resultado da Modelagem do Filtro Digital IIR bandstop.


Gerando os coeficientes para o Filtro IIR
 
Clicando no botão "Save C File"  o software (dsPIC FD Lite - Digital Filter Design) irá gerar um código fonte padrão na linguagem C para o filtro digital IIR construído, que pode ser visto no Código 1. Com ele você pode construir seu filtro DIgital IIR em outro compilador C para computador ou outro modelo de microcontrolador.


/***************************************************************************
****************************************************************************
*   File: E:\microchip_solutions_v2013-06-15\MyProjects\IIR FIlter\Digital Filter Design\Lowpass\lowpass1000Hz.c
*   Created by dsPIC FD Lite Version 1.0.0 Build 1 at 21:20:05 May 04 2015
*   C Code Generator - Version 4.0  
****************************************************************************
*  Code Fragment to implement filter
*
*  The functions defined in 'qed_filt.c' must be compiled and linked in.
*    This can be accomplished by either #include "qed_filt.c"
*    or by separately compiling and linking 'qed_filt.c'
*
*** following is actual code fragment
*  extern BiquadSections IIR_lowpass1000Hz;
*
*  init_biquad_float (&IIR_lowpass1000Hz);  // initialize filter structure 
*
*  IIR_lowpass1000Hz.filter ( x, y, n, &IIR_lowpass1000Hz);  // x is an array of input samples
*                                            // y is an array of output samples
*                                            // n is number of samples to process
*                                            // &IIR_lowpass1000Hz is a pointer to the
*                                            //    filter structure
*****************************************************************************
*  This is a complete program which can be compiled and run to test the filter.
*  To change this to a subroutine only, just add in this program or add globally
*     in "qed_cgen.h" the line with the definition of DEFINE_SUBROUTINE as follows
*  #define DEFINE_SUBROUTINE
*****************************************************************************
****************************************************************************/

/* qed_cgen.h contains definitions of filter structures and function prototypes */
#include "qed_cgen.h"

/* filter functions are in files 'qed_filt.c'  */

float lowpass1000Hz_num[  6] = {
   9.454345703125e-002F,  /* b[  1, 0] */
   1.891174316406e-001F,  /* b[  1, 1] */
   9.454345703125e-002F,  /* b[  1, 2] */
   6.185913085938e-002F,  /* b[  2, 0] */
   1.237487792969e-001F,  /* b[  2, 1] */
   6.185913085938e-002F}; /* b[  2, 2] */

float lowpass1000Hz_den[  6] = {
   1.000000000000e+000F,  /* a[  1, 0] */
  -8.160400390625e-001F,  /* a[  1, 1] */
   1.943054199219e-001F,  /* a[  1, 2] */
   5.000000000000e-001F,  /* a[  2, 0] */
  -5.340270996094e-001F,  /* a[  2, 1] */
   2.815856933594e-001F}; /* a[  2, 2] */

float lowpass1000Hz_m1[2];
float lowpass1000Hz_m2[2];

float lowpass1000Hz_gain = 1.000000000000e+000F; /* initial gain for cascade realization */
                          /* also applies to parallel realization */
float lowpass1000Hz_pars = 1.000000000000e+000F; /* scale value for parallel sections */

BiquadSections IIR_lowpass1000Hz = {

     2,     /* number of sections                     */
     0,     /* realization method                     */
            /*   0  Cascaded second order sections    */
            /*   1  Parallel second order sections    */
     1,     /* quantization: 0 off, 1 on              */
     1,     /* quantization type                      */
            /*   0  Floating point                    */
            /*   1  Fixed point fractional            */
     0,     /* realization type for cascaded sections only              */
            /*   0  Fixed point    - Transposed Canonic biquad sections */
            /*   1  Fixed point    - Canonic biquad sections            */
            /*   2  Floating point - 4 multiplies                       */
            /*   3  Floating point - 5 multiplies                       */
            /*   4  Floating point - recursive normal                   */
     0,     /* realization type for parallel sections only              */
            /*   0  Fixed point    - Transposed Canonic biquad sections */
            /*   1  Floating point - Transposed Canonic biquad sections */
   &lowpass1000Hz_gain,    /* pointer to gain for cascade/parallel realizations */
   &lowpass1000Hz_pars,    /* pointer to scale value for parallel sections      */
   lowpass1000Hz_num,      /* pointer to numerator coefficients                 */
   lowpass1000Hz_den,      /* pointer to denominator coefficients               */
   lowpass1000Hz_m1,       /* pointer to delay line 1                           */
   lowpass1000Hz_m2,       /* pointer to delay line 2                           */
   cas_blkfloat_fm1}; /* ptr to filter routine */

/* call the following function first and normally only once */
/* init_biquad_float (&IIR_lowpass1000Hz)  */
/*   where &IIR_lowpass1000Hz is a pointer to the BiquadSections */
/*   structure defining the filter */


/* call the following function to filter n samples */
/* IIR_lowpass1000Hz.filter (pIn, pOut, int n, &IIR_lowpass1000Hz); */

/*   where pIn  is a pointer to an array or buffer of samples to be filtered */
/*         pOut is a pointer to the array of filtered samples output by the filter */
/*         n    is the number of samples to filter */
/*         &IIR_lowpass1000Hz is a pointer to the structure defining the filter */


#ifndef DEFINE_SUBROUTINE

/* The following main program can be used to test the filter.         */
/*   input is in file 'in' and the filtered samples are in file 'out' */
/*   The input and output files are ascii floating point values       */
/*    e.g 1.0342 with 1 sample per line                               */
/* The input files can be created in DSPworks and exported as         */
/*   ascii floating point or any other system capable of creating     */
/*   ascii files with floating point values.                          */
/* The filtered output file can be imported into DSPworks as an ascii */
/*    floating point file and an FFT can be run to validate           */
/*    the frequency response.                                         */

#include "qed_filt.c"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>


#define INSZ1  1000
#define OUTSZ1 1000

static    float x[INSZ1], y[OUTSZ1];


int main(int argc, char *argv[])

{
  int i, in_count, file_status, error;

  FILE *fin;              /* input file of samples */
  FILE *fout;             /* output file of filtered samples */ 

  fprintf (stderr," ***** start filter test *****\n");

  fprintf (stderr," this program accepts 0,1 or 2 command line arguments\n");
  fprintf (stderr," the first  argument is the filename of the input file\n");
  fprintf (stderr," the second argument is the filename of the output file\n");
  fprintf (stderr," if there are 0 arguments, input and output is respectively\n");
  fprintf (stderr,"     stdin and stdout\n");
  fprintf (stderr," if only one argument is specified, then output is stdout\n");
  fprintf (stderr," if input is stdin rather than a file, then fscanf expects input\n");
  fprintf (stderr,"     from the console which may be piped in or entered directly\n");

  fin  = stdin;
  fout = stdout;
  error = 0;

  if (argc == 1) {
        fprintf(stderr," ***** waiting for input *****\n");
  }
  if (argc >= 2) {
      fin = fopen (argv[1], "r");
      if (fin == NULL) {
          fprintf(stderr,"\n error - Cannot open file %s for input\n", argv[1]);
          error = 1;
      }
  }
  if (argc >= 3) {
      fout = fopen (argv[2], "w");
      if (fout == NULL) {
          fprintf(stderr,"\n error - Cannot open file %s for output\n", argv[2]);
          error = 1;
      }
  }
  if (error) {
      fprintf(stderr," ***** end filter test *****\n");
      return(0);
  }
  
  

  init_biquad_float (&IIR_lowpass1000Hz);


  do {

      /* get input samples */ 
      for (in_count = 0; in_count < INSZ1; in_count++) { 
          file_status = fscanf(fin,"%f",&x[in_count]); 
          if (file_status != 1) 
              break; 
      }

      /* filter samples */ 

      if (in_count == 0) break;

      IIR_lowpass1000Hz.filter( x, y, in_count, &IIR_lowpass1000Hz);

      for (i = 0; i < in_count; i++)
          fprintf (fout,"%f\n",y[i]);

  } while (file_status == 1);

  fclose (fin); 
  fclose (fout);

  fprintf(stderr," ***** end filter test *****\n");
  return(1);

}
 
#endif 

Código 1 - Código C gerado com coeficientes para o Filtro IIR Lowpass.

O que nos interessa neste código são os valores dos coeficientes A1, A2, B1 e B2, que iremos utilizar no programa do PIC32MX. "A" e "B" são os pólos e zeros do Filtro.

Os valores dos coeficientes estão identificados em formato de matriz:
Coeficiente [ número da seção Biquadrada , número do coeficiente ]

Exemplo:
A [1,1] = coeficiente A1 da primeira seção biquadrada.
A [1,2] = coeficiente A2 da primeira seção biquadrada.
A [2,1] = coeficiente A1 da segunda seção biquadrada.
A [2,2] = coeficiente A2 da segunda seção biquadrada.

B [1,1] = coeficiente B1 da primeira seção biquadrada.
B [1,2] = coeficiente B2 da primeira seção biquadrada.
B [2,1] = coeficiente B1 da segunda seção biquadrada.
B [2,2] = coeficiente B2 da segunda seção biquadrada.

O IIR FIlter Biquad possui 2 pólos e 2 zeros para cada seção (seções de segunda ordem).

Para o Filtro IIR, A0 não é utilizado na função mips_iir16 ( ) da biblioteca DSP Library com o PIC32MX.

B0 não é necessário ser informado na função mips_iir16 ( ) da biblioteca DSP Library para o PIC32, pois no IIR FIlter Biquad, B0 = B2.

Por isto, no programa do PIC32MX, iremos utilizar somente os coeficientes A1, A2, B1 e B2 de cada seção.

Estamos modelando um Filtro Analógico Butterworth de 4ª ordem. Temos então 4 pólos e 4 zeros.

No menu Codegen / Microchip / dsPIC30 iremos gerar o código em Assembly com os coeficientes. Este código com os coeficientes para o filtro IIR lowpass é visto no Código 2.


; ..............................................................................
;    File   lowpass1000Hz.s
; ..............................................................................

        .equ lowpass1000HzNumSections, 2

; ..............................................................................
;
; Allocate and initialize filter coefficients
;
; These coefficients have been designed for use in the Transpose filter only

        .section .xdata

lowpass1000HzCoefs:
.hword    0x060D    ; b( 1,0)/2
.hword    0x0C1A    ; b( 1,1)/2
.hword    0x343A    ; a( 1,1)/2
.hword    0x060D    ; b( 1,2)/2
.hword    0xF391    ; a( 1,2)/2
.hword    0x07EB    ; b( 2,0)/2
.hword    0x0FD7    ; b( 2,1)/2
.hword    0x445B    ; a( 2,1)/2
.hword    0x07EB    ; b( 2,2)/2
.hword    0xDBF5    ; a( 2,2)/2

; ..............................................................................
; Allocate states buffers in (uninitialized) Y data space

        .section .ybss,  "b"

lowpass1000HzStates1:
        .space lowpass1000HzNumSections*2

lowpass1000HzStates2:
        .space lowpass1000HzNumSections*2

; ..............................................................................
; Allocate and intialize filter structure

        .section .data
        .global _lowpass1000HzFilter

_lowpass1000HzFilter:
.hword lowpass1000HzNumSections-1
.hword lowpass1000HzCoefs
.hword 0xFF00
.hword lowpass1000HzStates1
.hword lowpass1000HzStates2
.hword 0x0000

; ..............................................................................
; Sample assembly language calling program
;  The following declarations can be cut and pasted as needed into a program
;        .extern    _IIRTransposeFilterInit
;        .extern    _BlockIIRTransposeFilter
;        .extern    _lowpass1000HzFilter
;
;        .section    .bss
;
;     The input and output buffers can be made any desired size
;       the value 40 is just an example - however, one must ensure
;       that the output buffer is at least as long as the number of samples
;       to be filtered (parameter 4)
;input:        .space    40
;output:    .space    40
;        .text
;
;
;  This code can be copied and pasted as needed into a program
;
;
; Set up pointers to access input samples, filter taps, delay line and
; output samples.
;        mov    #_lowpass1000HzFilter, W0    ; Initalize W0 to filter structure
;        call    _IIRTransposeFilterInit    ; call this function once
;
; The next 4 instructions are required prior to each subroutine call
; to _BlockIIRTransposeFilter
;        mov    #_lowpass1000HzFilter, W0    ; Initalize W0 to filter structure
;        mov    #input, W1    ; Initalize W1 to input buffer 
;        mov    #output, W2    ; Initalize W2 to output buffer
;        mov    #20, W3    ; Initialize W3 with number of required output samples
;        call    _BlockIIRTransposeFilter    ; call as many times as needed


                           Código 2 - Código Assembly gerado com coeficientes para o Filtro IIR Lowpass.

A vantagem de utilizar os coeficientes do código em Assembly é o fato de já estar no formato fracionário. O programa do PIC32 necessita dos coeficientes no formato int16, então basta somente copiá-los para o nosso projeto. O código em C gera os coeficientes em float.

Este procedimento deve ser repetido para os 4 tipos de filtro IIR: lowpass, highpass, bandpass e bandstop. Os coeficientes A1, A2, B1 e B2 são diferentes para cada design  elaborado.


O Hardware

Neste projeto utilizaremos:
a) PIC32MX USB Starter Kit II;
b) I/O Expansion Board;
c) GFXv3 Graphic Display.

O PIC32 irá gerar um sinal senoidal de frequência variável, aplicar o filtro digital IIR, calcular a Transformada Rápida de Fourier, a fim de verificar a resposta de cada tipo de filtro e comparar com os dados teóricos obtidos com o dsPIC FD Lite e mostrar no display os respectivos gráficos.



O código fonte do PIC32

 No terceiro passo elaboraremos o programa do PIC no compilador MPLAB C32 da Microchip com o ambiente de desenvolvimento MPLAB X. O Programa foi desenvolvido em linguagem C.

Incluímos em nosso projeto a biblioteca "dsplib_dsp.h" que possui as funções para o cálculo dos filtros digitais IIR, FIR e outras utilizadas em cálculos de processamentos digitais, como a FFT.

O Código 3 mostra a inclusão da biblioteca e a declaração dos vetores do Filtro Digital IIR.




#include "dsplib_dsp.h"   // library to FFT and IIR FIlter functions


// vetores para iir filter
#define B  2    //2 Biquad Sections
biquad16 bq[B];
int16 coeffs[4*B];
int16 delayline[2*B];
int16 indata, outdata;
int16 sinal_Filtrado [N*2];


Código 3 - Incluindo a biblioteca DSP e declarando os vetores.

Após o início do programa, na função main ( ) temos que iniciar as variáveis e atribuir os coeficientes do filtro. Os valores para A1, A2, B1 e B2 foram retirados do código Assembly gerado pelo dsPIC Filter Design Lite.

A inicialização das variáveis e coeficientes são vistos no Código 4.

Os coeficientes são carregados para A1, A2, B1 e B2 de acordo com o tipo de filtro selecionado: lowpass, highpass, bandpass ou bandstop.




int main(void)
{
   ...


    // inicializa componentes de hardware e periféricos
    InitializeBoard();

    while(1)
    {

        // carrega coeficientes de acordo com o filtro selecionado
        iir_coefficients_setup ();

    ...
    ...
    ...
    }
}

void iir_coefficients_setup ()
{
    int i;
        //inicializa com zero delayline do iir filter
        for (i = 0; i < 2*B; i++)
            delayline[i] = 0;
   
    
        /* altera coeficientes do filtro IIR conforme a selecao */
        // load coefficients into bq here
        switch  (filterType) {
            case FILTER_OFF:
                
                break;
                
            case LOWPASS:
                bq [0].a1 = 0x343A; //a( 1,1)/2
                bq [0].a2 = 0xF391; //a( 1,2)/2
                bq [0].b1 = 0x0C1A; //b( 1,1)/2
                bq [0].b2 = 0x060D; //b( 1,2)/2

                bq [1].a1 = 0x445B; //a( 2,1)/2
                bq [1].a2 = 0xDBF5; //a( 2,2)/2
                bq [1].b1 = 0x0FD7; //b( 2,1)/2
                bq [1].b2 = 0x07EB; //b( 2,2)/2
                break;
                
            case HIGHPASS:
                bq [0].a1 = 0xFBAF; //a( 1,1)/2
                bq [0].a2 = 0xFD67; //a( 1,2)/2
                bq [0].b1 = 0xE0DC; //b( 1,1)/2
                bq [0].b2 = 0x0F92; //b( 1,2)/2

                bq [1].a1 = 0xFA00; //a( 2,1)/2
                bq [1].a2 = 0xE360; //a( 2,2)/2
                bq [1].b1 = 0xD4B0; //b( 2,1)/2
                bq [1].b2 = 0x15A7; //b( 2,2)/2
                break;
             case BANDPASS:
                bq [0].a1 = 0xFF3C; //a( 1,1)/2
                bq [0].a2 = 0xD061; //a( 1,2)/2
                bq [0].b1 = 0x0000; //b( 1,1)/2
                bq [0].b2 = 0xF7D0; //b( 1,2)/2

                bq [1].a1 = 0x1FDD; //a( 2,1)/2
                bq [1].a2 = 0xCFC9; //a( 2,2)/2
                bq [1].b1 = 0x0000; //b( 2,1)/2
                bq [1].b2 = 0xEFDC; //b( 2,2)/2
                break;
                
            case BANDSTOP:
                bq [0].a1 = 0x91BB; //a( 1,1)/2
                bq [0].a2 = 0xCF87; //a( 1,2)/2
                bq [0].b1 = 0xFFD1; //b( 1,1)/2
                bq [0].b2 = 0x0102; //b( 1,2)/2

                bq [1].a1 = 0x7135; //a( 2,1)/2
                bq [1].a2 = 0xCD3B; //a( 2,2)/2
                bq [1].b1 = 0xECE9; //b( 2,1)/2
                bq [1].b2 = 0x6872; //b( 2,2)/2
                break;  
                
            default :
                break;              
       }               


        mips_iir16_setup(coeffs, bq, B); 
    
}

Còdigo 4 - Carregando os coeficientes e inicializando a estrutura do Filtro IIR.

Primeiro deve-se iniciar com zero o vetor delayline [ ].

Depois carregam-se os valores dos coeficientesem bq [ ], de acordo com o tipo do filtro IIR selecionado.

Agora então podemos chamar a função mips_iir16_setup ( ) da biblioteca DSP do PIC32MX que se encarrega de inicializar a estrutura do filtro IIR.

O Código 5 mostra a aplicação do filtro propriamente dito em uma amostra do sinal de entrada.



/*  calculate iir filter biquad */
void iir_filter ()
{
    int i;

    for (i=0;i<N*2;i++)
    {
       // get input data value into indata
        indata = sine_wave_input[i] ;

        //calculate iir filter biquad
        outdata = mips_iir16(indata, coeffs, delayline, B, 1);

        sinal_Filtrado [i] = outdata;
    } //end for
}

Código 5 - Função que aplica o FIltro IIR ao sinal de entrada.

A função responsável pelo cálculo do Filtro IIR Biquad é a mips_iir16 ( ), e pertence à DSP Library da Microchip para o PIC32MX.

Nosso programa executa o Código 6 em um loop infinito: gera uma onda senoidal de frequência variável; aplica o filtro IIR ao sinal gerado; calcula a FFT ao sinal Filtrado e atualiza o display com os gráficos.


            sineWaveGeneration (signal_freq, 7680.0 );   //gera sinal senoidal
            iir_filter ();              //aplica o filtro iir ao sinal gerado
            moveVectorInputFFT ();      //move sinal gerado para entrada do calculo da FFT
            fftCalculate ();            //calcula a FFT do sinal gerado
            updateDisplay ();           //atualiza a tela


Código 6 - Loop do programa no código main.

O resultado final do nosso projeto é visto na Figura 14. Interligando um display gráfico QVGA 320x240 vemos o sinal gerado sendo aplicado ao Filtro Digital IIR.

Figura 14 - Sinais sendo mostrados no display gráfico.


A Figura 15 mostra a resposta na prática para o nosso filtro passa-baixa de 1000 Hz.

Figura 15 - Resposta do IIR FIlter lowpass.

A Figura 16 mostra a resposta na prática para o nosso filtro passa-alta de 2000 Hz.

Figura 16 - Resposta do IIR FIlter highpass.

A Figura 17 mostra a resposta na prática para o nosso filtro passa-banda de 1500 Hz à 2000 Hz.

Figura 17 - Resposta do IIR FIlter bandpass.

A Figura 18 mostra a resposta na prática para o nosso filtro rejeita-banda de 2000 Hz à 2500 Hz

Figura 18 - Resposta do IIR FIlter bandstop.

Estas curvas podem ser comparadas com aquelas do início deste post, geradas pelo dsPIC FD Lite, a fim de verificar a performance do Filtro Digital IIR da DSP Library com o PIC32MX!

Lembrando que a frequência de amostragem é 7680 Hz e o sinal gerado varia de 0 à 3840 Hz, respeitrando o Teorema de Nyquist.

O Vídeo 1 mostra o projeto do Filtro Digital IIR elaborado no MPLAB X  e rodando com o PIC32MX.

Vídeo 1 - IIR FIlter com PIC32MX.


Os cálculos da Fast Fourier Transform (FFT) para o PIC32MX você encontra aqui.




Referências:

ORFANIDS, Sophocles J. Introduction to Signal Processing. Pearson Education/Prentice Hall/Sophocles J. Orfanidis.2009. Disponível em: http://www.ece.rutgers.edu/~orfanidi/intro2sp

PALACHERLA, Amar. Implementing IIR Digital Filters. Application Note AN540. Microchip Technology Inc.1997.

RAMU, B. K. Anantha. Implementing FIR and IIR Digital Filters using PIC18 microcontrollers. Application Note AN852. Microchip Technology Designs. India. 2002.

SKANDA, Vinaya. Sine Wave Generator Using Numerically Controller PIC. Application Note AN1523A. Microchip Technology Inc. 2013.

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

MANCHESTER. Z-transforms & IIR-type digital filters.  University of Manchester. School of Computer Science. Digital Media Processing. Section 5. 2009. 19p. Disponível em: www.cs.man.ac.uk/~barry/mydocs/COMP30291/Notes/DMP09-5.doc

MICROCHIP.  Using the Fast Fourier Transform (FFT) for Frequency Detection. Code Examples CE018. Microchip Technology Inc. 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



Anexo A

Código básico para aplicar o cálculo de 1 tipo de Filtro Digital IIR.



#include "dsplib_dsp.h"   // library to FFT and IIR FIlter functions


// vetores para iir filter
#define B  2    //2 Biquad Sections
biquad16 bq[B];
int16 coeffs[4*B];
int16 delayline[2*B];
int16 indata, outdata;
int16 sinal_Filtrado [N*2];



int main(void)
{
   ...


    // inicializa componentes de hardware e periféricos
    InitializeBoard();

    
    // carrega coeficientes de acordo com o filtro selecionado
    iir_coefficients_setup ();


    while(1)
        {

        iir_filter ();    //aplica o filtro iir ao sinal gerado

        ...
        ...
        ...
        }
}

void iir_coefficients_setup ()
{
    int i;
        //inicializa com zero delayline do iir filter
        for (i = 0; i < 2*B; i++)
            delayline[i] = 0;
   
    
        // load coefficients into bq here
        
                bq [0].a1 = 0x343A; //a( 1,1)/2
                bq [0].a2 = 0xF391; //a( 1,2)/2
                bq [0].b1 = 0x0C1A; //b( 1,1)/2
                bq [0].b2 = 0x060D; //b( 1,2)/2

                bq [1].a1 = 0x445B; //a( 2,1)/2
                bq [1].a2 = 0xDBF5; //a( 2,2)/2
                bq [1].b1 = 0x0FD7; //b( 2,1)/2
                bq [1].b2 = 0x07EB; //b( 2,2)/2
 
        mips_iir16_setup(coeffs, bq, B); 
    
}



/*  calculate iir filter biquad */
void iir_filter ()
{
    int i;

    for (i=0;i<N*2;i++)
    {
       // get input data value into indata
        indata = sine_wave_input[i] ;

        //calculate iir filter biquad
        outdata = mips_iir16 (indata, coeffs, delayline, B, 1);

        sinal_Filtrado [i] = outdata;
    } //end for
}





domingo, 15 de março de 2015

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



sábado, 7 de março de 2015

Como Utilizar a FFT com o dsPIC - DSP Library

Este post mostra um exemplo de como utilizar a Fast Fourier Transformer (FFT) da Microchip DSP Library da família dsPIC.


dsPICs  são Controladores Digitais de Sinais de 16-bit da Microchip, uma junção de Microcontrolador com Processador Digital de Sinais (DSP). Possui a memória RAM dividida em 2 partes: Memória X para acesso das instruções da parte do Microcontrolador e Memória Y para acesso das instruções de DSP. Os cálculos de DSP são executados com variáveis do tipo fracionária (fractional Q15), com range entre -1.0 e +1.0. 

As instruções de DSP são tipicamente utilizadas para o Processamento Digital de Sinais. Como todo DSP, o dsPIC também possui a instrução MAC!

A DSP Library da Microchip possui diversas funções como:
- operação com vetores: máximo, mínimo,  soma, subtração, multiplicação, convolução, correlação, etc.;
- operação com janelas: Bartlett, Blackman, Hamming, Kaiser, etc;
- operação com matrizes;
- filtros digitais FIR e IIR;
- transformadas FFT e IFFT;

Aqui iremos apresentar a utilização da função FFTComplexIP ( ) da DSP Library utilizando a linguagem C. O cálculo da FFT da DSP Library está escrito em Assembly, pois utiliza as instruções de DSP, para permitir o menor tempo de processamento.


Sinal de Entrada para o Cálculo

Para o sinal de entrada que será aplicado a FFT iremos gerar uma onda quadrada com o auxílio do Software dsPICworks, da Microchip.

No menu Generator encontramos a opção Square... para gerar uma onda quadrada.

A Figura 1 mostra os parâmetros que devemos preencher para o software gerar 128 amostras de uma onda quadrada com a frequência de 780 Hz, taxa de amostragem de 7680 sps, do tipo fracionário de 16 bit.

FIgura 1 - Configurações para gerar uma onda quadrada de 780 Hz e 7680 sps no dsPICworks.


A onda quadrada gerada, no domínio do tempo, é vista na Figura 2.

Figura 2 - Onda gerada no dsPICWork.


No menu File, opção Export File... podemos exportar os dados da onda gerada pelo dsPICworks.

O vetor com os dados gerados já adicionado em nosso programa está contido no Código 1. Repare que o vetor é do tipo fractional, 16-bit, com 128 dados de valores -1.0 e +1.0.



fractional sinal_Entrada [128]= {            //vetor que ira armazenar os dados de entrada
                                          //valores iniciais declarados para simulacao
                                          //e teste do programa
                                          //square signal 780 Hz 7.68 kSPS
                                          //sinal gerado com o dsPICWorks
    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, 0x8001, 0x8001,
    0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
    0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
    0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
    0x8001, 0x8001, 0x8001};

Código 1 - Vetor de dados da onda quadrada exportada pelo dsPICworks para ser utilizado no MPLAB.


Dentro do arquivo header DSP.h da Biblioteca DSP podem ser encontrados detalhes dos argumentos da função FFTComplexIP ( ). Quem tiver interesse pode encontrá-lo no Código 2.


/****************************************************************************
*
* Interface to transform operations.
*
* A set of linear discrete signal transformations (and some of the inverse
* transforms) are prototyped below. The first set applies a Discrete Fourier
* transform (or its inverse) to a complex data set. The second set applies
* a Type II Discrete Cosine Transform (DCT) to a real valued sequence.
*
* A complex valued sequence is represented by a vector in which every pair
* of values corresponds with a sequence element. The first value in the pair
* is the real part of the element, and the second its imaginary part (see
* the declaration of the 'fractcomplex' structure at the beginning of this
* file for further details). Both, the real and imaginary parts, are stored
* in memory using one word (two bytes) each, and must be interpreted as Q.15
* fractionals.
*
* The following transforms have been designed to either operate out-of-place,
* or in-place. The former type populates an output sequence with the results
* of the transformation. In the latter, the input sequence is (physically)
* replaced by the transformed sequence. For out-of-place operations, the user
* must provide with enough memory to accept the results of the computation.
* The input and output sequences to the FFT family of transforms must be
* allocated in Y-Data memopry.
*
* The transforms here described make use of transform factors which must be
* supplied to the transforming function during its invokation. These factors,
* which are complex data sets, are computed in floating point arithmetic,
* and then transformed into fractionals for use by the operations. To avoid
* excessive overhead when applying a transformation, and since for a given
* transform size the values of the factors are fixed, a particular set of
* transform factors could be generated once and used many times during the
* execution of the program. Thus, it is advisable to store the factors
* returned by any of the initialization operations in a permanent (static)
* vector. The factors to a transform may be allocated either in X-Data or
* program memory.
*
* Additional remarks.
*
* A) Operations which return a destination vector can be nested, so that
*    for instance if:
*
*    a = Op1 (b, c), with b = Op2 (d), and c = Op3 (e, f), then
*
*    a = Op1 (Op2 (d), Op3 (e, f))
*
****************************************************************************/

/* Transform operation prototypes. */


extern fractcomplex* TwidFactorInit (    /* Initialize twiddle factors */
                    /* WN(k) = exp(i*2*pi*k/N) */
                    /* computed in floating point */
                    /* converted to fractionals */
   int log2N,                /* log2(N), N complex factors */
                       /* (although only N/2 are computed */
                       /* since only half of twiddle factors */
                    /* are used for I/FFT computation) */
   fractcomplex* twidFactors,        /* ptr to twiddle factors */
   int conjFlag                /* indicates whether to generate */
                    /* complex conjugates of twiddles */
                    /* 0 : no conjugates (default) */
                    /* 1 : conjugates */

                    /* twidfact returned */
                    /* only the first half: */
                    /* WN(0)...WN(N/2-1) */
                    /* (or their conjugates) */
);

/*...........................................................................*/

extern fractcomplex* BitReverseComplex (    /* Bit Reverse Ordering */
                    /* (complex) */
   int log2N,                /* log2(N), N is vector length */
   fractcomplex* srcCV            /* ptr to source complex vector */
                    /* MUST be N modulo aligned */

                    /* srcCV returned */
);

/*...........................................................................*/

extern fractcomplex* FFTComplex (    /* Fast Fourier Transform */
                    /* (complex, out-of-place) */
   int log2N,                /* log2(N), N-point transform */
   fractcomplex* dstCV,            /* ptr to destination complex vector */
                       /* with time samples */
                    /* in natural order */
                    /* MUST be N modulo aligned */
   fractcomplex* srcCV,            /* ptr to source complex vector */
                       /* with time samples */
                    /* in natural order */
   fractcomplex* twidFactors,        /* base address of twiddle factors */
                    /* either in X data or program memory */
                    /* if in X data memory, it points at */
                       /* WN(0).real */
                    /* if in program memory, base is the */
                    /* offset from program page boundary */
                    /* to address where factors located */
                    /* (inline assembly psvoffset ()) */
   int factPage                /* if in X data memory, set to */
                       /* defined value COEFFS_IN_DATA */
                    /* if in program memory, page number */
                    /* where factors are located */
                    /* (inline assembly psvpage ()) */

                    /* dstCV returned */
                       /* with frequency components */
                    /* in natural order */
                    /* and scaled by 1/(1<<log2N) */
);

/*...........................................................................*/

extern fractcomplex* FFTComplexIP (    /* Fast Fourier Transform */
                    /* (complex, in-place) */
   int log2N,                /* log2(N), N-point transform */
   fractcomplex* srcCV,            /* ptr to source complex vector */
                       /* with time samples */
                    /* in natural order */
   fractcomplex* twidFactors,        /* base address of twiddle factors */
                    /* either in X data or program memory */
                    /* if in X data memory, it points at */
                       /* WN(0).real */
                    /* if in program memory, base is the */
                    /* offset from program page boundary */
                    /* to address where factors located */
                    /* (inline assembly psvoffset ()) */
   int factPage                /* if in X data memory, set to */
                       /* defined value COEFFS_IN_DATA */
                    /* if in program memory, page number */
                    /* where factors are located */
                    /* (inline assembly psvpage ()) */

                    /* srcCV returned */
                       /* with frequency components */
                    /* in bit reverse order */
                    /* and scaled by 1/(1<<log2N) */
);
Código 2 - Parte do arquivo DSP.h da Microchip DSP Library com declaração das funções da FFT.


No início do nosso projeto devemos incluir o header DSP.h da biblioteca DSP, conforme mostra o Código 3.

O arquivo Header do dsPIC30F4013 também é incluso aqui, aliás, este é o DSC que iremos utilizar para a demonstração da FFT. Possui capacidade para 30 MIPS @117 MHz, 48 kB de Flash para memória de programa, 1 kB de memória RAM X e 1 kB de memória RAM Y, 1 kB de EEPROM e 13 canais ADC 12-bit com 200 ksps.

#include "dsp.h"
#include "p30f4013.h"
#include <stdio.h> 
Código 3 - Incluindo a biblioteca DSP.h ao projeto.

Os seguintes arquivos de biblioteca (.a) foram adicionados ao projeto:
- libc-coff.a;
- libdsp-coff.a;
- libm-coff.a;
- libp30f4013.a;
- libpic30-coff.a.


A Transformada Rápida de Fourier

Bom, agora temos que compreender alguns conceitos básicos da FFT. Para mais detalhes você pode acessar este outro post: http://mrgptu.blogspot.com.br/2015/01/analisador-de-espectro-fft-pic.html

A Transformada de Fourier (FT) é utilizada para converter uma função ou sinal do domínio do tempo para o domínio da frequência e vice-versa.

A Transformada Discreta de Fourier (DFT) é aplicada para sinais discretos no tempo (dados coletados por um conversor analógico digital por exemplo).

A Transformada Rápida de Fourier (FFT) foi desenvolvida para uma execução mais rápida dos cálculos da FT possibilitando o seu uso na computação. Dentre as suas principais aplicações estão a filtragem digital, reconhecimento de padrões, bargraph para sinais de áudio, equalizadores, eliminações de ruído e interferência de imagens.

A FFT representa a soma de uma série de ondas senoidais de diferentes frequências, fases e amplitudes.

Para entender a diferença do sinal analisado no domínio do tempo e no domínio da frequência vamos a um exemplo prático utilizando o software dsPICworks, da Microchip. Com ele geramos alguns sinais. 

Na Figura 3a temos um sinal senoidal de 60 Hz (domínio do tempo, eixo x do gráfico).
 

  Figura 3a - Sinal Senoidal de 60 Hz no tempo.


 Na Figura 3b um sinal senoidal de 600 Hz (domínio do tempo).


 Figura 3b - Sinal Senoidal de 600 Hz no tempo.


Efetuando a soma dos sinais senoidais A e B obtemos o sinal C mostrado na Figura 3c (domínio do tempo).



  Figura 3c - Soma dos Sinais Senoidais de 60 Hz e 600 Hz no tempo.


O resultado da Transformada Rápida de Fourier aplicada ao sinal C é exibido na Figura 3D. O sinal agora é visto no domínio da frequência (eixo x do gráfico). Observe 2 picos no gráfico, que representam a energia do espectro de 60 Hz e de 600 Hz.

 Figura 1d - Resultado da FFT aplicada ao sinal.


A análise de um sinal no domínio da frequência se torna muito interessante e/ou mais amigável conforme a aplicação, por exemplo, em um VU Meter, equalizador gráfico de áudio ou na análise das frequências harmônicas da rede de energia elétrica (múltiplas de 60 Hz).


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


//---------------------- Define constantes -------------------------------
// constantes da FFT
#define PONTOS_FFT    128         //numero de pontos da FFT
#define LOG2_ESTAGIOS_BUTTERFLY 7    //numero de estagios "Butterfly" da FFT
#define TAXA_AMOSTRAGEM        7680    //taxa de amostragem do sinal de entrada 
Código 4 - 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,680 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).


O Código Fonte

Retornando ao nosso projeto, iremos declarar as outras variáveis e vetores de dados que serão utilizados, Veja no Código 5 que o vetor sinalComplexo é do tipo fractcomplex (fracionário com partes real e parte imaginária). Este vetor irá conter os dados de entrada complexos para o cálculo da FFT. Este vetor é alocado na memória RAM Y, acessada pelas instruções de DSP contidas na função FFT.




int    picoFrequenciaBin = 0;                
unsigned long picoFrequencia = 0;

fractional *temp3_pnt, tempFrac;

// Entrada do sinal para a FFT declarada na memoria Y
fractcomplex sinalComplexo[PONTOS_FFT]__attribute__ ((section (".ydata, data, ymemory"),aligned (PONTOS_FFT * 2 *2))); 

fractional SquareMagnitude [PONTOS_FFT/2];
Código 5 - Declara variáveis e vetores.


Os cálculos da FFT utilizam coeficientes chamados Twiddle. Estes coeficientes podem ser calculados através da equação:

                                 WN (kn) = exp [-(j*2*pi*k*n)/N]

No nosso programa utilizamos os coeficientes declarados na memória de programa. A quantidade de coeficientes depende do número de pontos da FFT. O Código 6 contém os coeficientes Twiddle para 128 pontos da FFT. Os coeficientes também podem estar alocados na memória RAM, o que permite o processamento da FFT em menor tempo, mas nem sempre temos espaço em memória RAM disponível (1 kB para o dsPIC30F4013), então podemos optar por declarar os coeficientes na memória de programa (48 kB para o dsPIC30F4013).

// ------------- Declara constantes gravadas na flash -------------------------

// declara vetor dos coeficientes Twiddle da FFT na memoria de programa
// Define 128 constantes "twiddle" da FFT: WN (kn) = exp [-(j*2*pi*k*n)/N]
const fractcomplex twiddleFactors[] __attribute__ ((space(prog), aligned (PONTOS_FFT*2)))=
        {
        0x7FFF, 0x0000, 0x7FD9, 0xF9B8, 0x7F62, 0xF374, 0x7E9D, 0xED38,
        0x7D8A, 0xE707, 0x7C2A, 0xE0E6, 0x7A7D, 0xDAD8, 0x7885, 0xD4E1,
        0x7642, 0xCF04, 0x73B6, 0xC946, 0x70E3, 0xC3A9, 0x6DCA, 0xBE32,
        0x6A6E, 0xB8E3, 0x66D0, 0xB3C0, 0x62F2, 0xAECC, 0x5ED7, 0xAA0A,
        0x5A82, 0xA57E, 0x55F6, 0xA129, 0x5134, 0x9D0E, 0x4C40, 0x9930,
        0x471D, 0x9592, 0x41CE, 0x9236, 0x3C57, 0x8F1D, 0x36BA, 0x8C4A,
        0x30FC, 0x89BE, 0x2B1F, 0x877B, 0x2528, 0x8583, 0x1F1A, 0x83D6,
        0x18F9, 0x8276, 0x12C8, 0x8163, 0x0C8C, 0x809E, 0x0648, 0x8027,
        0x0000, 0x8000, 0xF9B8, 0x8027, 0xF374, 0x809E, 0xED38, 0x8163,
        0xE707, 0x8276, 0xE0E6, 0x83D6, 0xDAD8, 0x8583, 0xD4E1, 0x877C,
        0xCF04, 0x89BE, 0xC946, 0x8C4A, 0xC3A9, 0x8F1D, 0xBE32, 0x9236,
        0xB8E3, 0x9592, 0xB3C0, 0x9931, 0xAECC, 0x9D0E, 0xAA0A, 0xA129,
        0xA57E, 0xA57E, 0xA129, 0xAA0A, 0x9D0E, 0xAECC, 0x9931, 0xB3C0,
        0x9592, 0xB8E3, 0x9236, 0xBE32, 0x8F1D, 0xC3A9, 0x8C4A, 0xC946,
        0x89BE, 0xCF04, 0x877C, 0xD4E1, 0x8583, 0xDAD8, 0x83D6, 0xE0E6,
        0x8276, 0xE707, 0x8163, 0xED38, 0x809E, 0xF374, 0x8027, 0xF9B8
        } ;

Código 6 - Declara coeficientes Twiddle na memória de programa.


O Código 7 mostra o início da função main. De imediato criamos e declaramos algumas variáveis do tipo ponteiro.

int main (void)            //inicio do programa
{
    int i = 0;
    fractional *p_real = &sinalComplexo[0].real ;
    fractcomplex *p_complexo = &sinalComplexo[0] ;
    fractional *p_fract = &sinalComplexo[0].real ;

    ...

Código 7 - Início do programa main e declaração de ponteiros para os vetores de dados.


No Código 8 inicia-se o loop while. Encontramos 3 loops for:
1) movemos o vetor de entrada (a onda quadrada gerada pelo dsPICworks para teste do cálculo da FFT) para o vetor complexo utilizado pela FFT;
2) escalonamos o sinal de entrada de -1,0 a +1,0 para -0,5 a +0,5, para evitar overflow dos cálculos;
3) organizamos as partes real e imaginária dos dados do vetor complexo utilizado pela FFT. A parte real contém os dados do sinal gerado, ou de um ADC. A parte imaginária neste caso é zero.


  while (1)      //loop infinito
  {

     
    p_real = &sinalComplexo[0].real ;              //inicializa ponteiro
    p_complexo = &sinalComplexo[0];                //inicializa ponteiro
    
    /* Move buffer do sinal de entrada (real)para vetor da FFT    */
    for ( i = 0; i < PONTOS_FFT; i++ )             
    {
        *p_real = sinal_Entrada [i] ;        
        *p_real++;                                 //incrementa o ponteiro     
    }

    
    /* Escalona o vetor para o range [-0.5, +0.5]                */
    p_real = &sinalComplexo[0].real ;             //inicializa ponteiro
    for ( i = 0; i < PONTOS_FFT; i++ )         
    {
        *p_real = *p_real >>1 ;                     //desloca 1 bit para a direita
        *p_real++;            
    }                    


    /* Converte vetor real para vetor complexo                    */
    p_real = &sinalComplexo[(PONTOS_FFT/2)-1].real ;    //inicializa o ponteiro para parte real do vetor  
    p_complexo = &sinalComplexo[PONTOS_FFT-1] ; 

    for ( i = PONTOS_FFT; i > 0; i-- )         
    {
        (*p_complexo).real = (*p_real--);    
        (*p_complexo--).imag = 0x0000;            //Não possui parte imaginaria, escreve valor zero
    }

    ...

Código 8 - Prepara vetor de dados complexos de entrada para cálculos da FFT.


Agora podemos aplicar a Transformada Rápida de Fourier, a FFT, na onda quadrada de 780 Hz!

A função FFTComplexIP se encarrega desta tarefa!

Você encontra a chamada da FFT no Código 9. Após o cálculo da FFT é necessário organizar os dados de saída com a função BitReverseComplex.

    /* Cálcula a FFT - Transformada Rapida de Fourier             */
    FFTComplexIP (LOG2_ESTAGIOS_BUTTERFLY, &sinalComplexo[0], (fractcomplex *) __builtin_psvoffset(&twiddleFactors[0]), (int) __builtin_psvpage(&twiddleFactors[0]));


    /* Organiza dados da saída da FFT - Bit-reverso                */
    BitReverseComplex (LOG2_ESTAGIOS_BUTTERFLY, &sinalComplexo[0]);

Código 9 - Cálculo da FFT e organização do vetor de dados de saía.


No Código 10 você encontra um programa adicional para calcular a soma do quadrado das partes real e imaginária, encontrar a frequência com maior espectro e o cálculo do seu módulo. Isto é útil para criar um gráfico para mostrar o resultado da FFT.

    /* Calcula a soma do quadrado das partes real e imaginária 
       transformando vetor complexo em um vetor real 
            Sinal = real^2 + imaginario^2                      */
    SquareMagnitudeCplx(PONTOS_FFT, &sinalComplexo[0], &sinalComplexo[0].real);

    /* Move buffer do vetor complexo para vetor comum             */
    p_fract = &sinalComplexo[0].real ;                //inicia ponteiro
    for ( i = 0; i < PONTOS_FFT/2; i++ )             
    {
        SquareMagnitude [i] = *p_fract;        
        *p_fract++;                                     //incrementa o ponteiro                                               
    }


    /* Encontra a frequency Bin ( indice da frequencia) que possui a maior energia */
    VectorMax(PONTOS_FFT/2, &SquareMagnitude[0], &picoFrequenciaBin);

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

    /* Encontra  o modulo da frequencia de pico                    */
    temp3_pnt = &SquareMagnitude[0] ;      //inicializa ponteiro    
    temp3_pnt = temp3_pnt + picoFrequenciaBin;
    tempFrac = *temp3_pnt;
    tempFloat = tempFrac;
    tempFloat = tempFloat / 32767;            // Converte fracionario para float
    modulo = sqrt (tempFloat);                // Extrai a raiz quadrada da soma dos quadrados da parte real/imaginaria

Código 10 - Calcula módulo e encontra frequência de maior energia.


A Simulação dos cálculos no MPLAB IDE

A Figura 4 mostra o projeto aberto no MPLAB IDE v8.83.

Figura 4 - Projeto no MPLAB IDE v8.83 com compilador MPLAB C30.


Na Figura 5 é visualizada a janela Watch da simulação no MPLAB contendo as variáveis e vetores do programa após serem realizados os cálculos.


Figura 5 - Tabela Watch com resultado do cálculo da FFT sobre o sinal quadrado gerado.


Observe que a variável picoFrequencia mostra o valor de 780 [Hz]. Esta é a frequência da nossa onda quadrada gerada no dsPICwork, lembra?

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 código completo deste programa você encontra no Anexo A.

O Vídeo 1 mostra este projeto rodando no MPLAB IDE com o dsPIC30F4013. É simulado passo-a-passo a execução de cada função, incluindo a FFTComplexIP que efetua o cálculo da Fast Fourier Transform - FFT,  e as alterações das variáveis na tabela Watch.

Vídeo 1 - Exemplo da FFT rodando no MPLAB IDE


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

Vídeo 2 - VU Meter Audio Spectrum Analyser com o dsPIC30F.



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

A FFT é testada com um gerador de sinal aplicando uma onda senoidal variando de 20 Hz à 1,5 kHz.

Para desenvolver sua própria aplicação e não utilizar mais os dados gerados pelo dsPICwork, você deve configurar o conversor analógico-digital para amostrar o sinal em uma taxa de 7680 sps (7680 Hz) e no formato fracionário com sinal (ADC_FORMAT_SIGN_FRACT). 128 amostras do ADC serão copiadas para o vetor sinal_entrada [ ]. Pronto, agorá é só executar o código aqui sugerido.

Outros valores de Fs, pontos da FFT e Fo podem ser utilizados. Quanto maior Fs (frequencia de amostragem) mais rápido terá que ser seu ADC. Quanto mais pontos da FFT mais memória RAM é necessária para armazenar os vetores de dados e maior será o processamento da FFT (quantidade maior de cálculos).

Um exemplo para a aplicação da FFT você encontra neste post: http://mrgptu.blogspot.com.br/2015/01/analisador-de-espectro-fft-pic.html

Trata-se de um Analisador de Espectro que analisa as harmônicas da rede de energia elétrica (60 Hz). As funções da FFT são as mesmas utilizadas neste artigo, inclusive mesma Fs, pontos da FFT e Fo.


Anexo A - Código completo da FFT com o dsPIC30F4013

/************************************************************************************************
********************* DEMONSTRAÇÃO DO USO DA FFT NO dsPIC ***************************************
*************************************************************************************************
  
Caracteristicas:
- Microcontrolador & DSP dsPIC30F4013 trabalhando a 30 MIPS, com oscilador interno, WDT desligado
- Analise do espectro de frequencia atraves da Transformada Rapida de Fourier - FFT
- Calcula a frequencia com maior espectro e seu modulo 
- Sinal gerado para simulacao dos calculos: square signal, 780 Hz, 7680 SPS, amplitude -1.0 a +1.0

- Dados da FFT:
    Pontos: 128; 128 amostragens no tempo transformarao em 128 pontos de frequencia
     fs: 7,68 kHz; taxa de amostragem
     fo: 60 Hz; passo entre as frequencias       



Rev. 1.2 Criada aplicação de exemplo para uso da FFT


;##### M.R.G. ######################################
;##### MARÇO/2015 ##################################
;##### USADO MPLAB IDE V8.83 E MPLAB C30 v2.02 #####

;*********** CLOCK INTERNO DE 7.37 MHZ x PLL16 30 MIPS *******************************************

;***** DATA SPACE MEMORY MAP 2K MSB    LSB *******************************************************
;    X DATA RAM :             0X8001 0X0800           A        0X0BFF 0X0BFE instruçóes MCU e DSC
;    Y DATA RAM :             0X0C01 0X0C00           A        0X0FFF 0X0FFE instruções DSC
;***************** obs.: 0x0800 contém endereço do TOS (topo da pilha) ***************************
;*************************************************************************************************  
;***** EEPROM. 1 Kbytes of data EEPROM with an address globalrange from 0x7FFC00 to 0x7FFFFE *****
;*************************************************************************************************    
*/
 
//------------------ Inclui ficheiros auxiliares -------------------------
#include "dsp.h"
#include "p30f4013.h"
#include <stdio.h> 


/*-------------------------- FUSES ------------------------------*/
//--------- BITS DE CONFIGURACAO DO DISPOSITIVO -----------------
_FOSC  (0xBFE3); //(CSW_FSCM_OFF & FRC_PLL16);         //30 MIPS
//_FOSC  (0xB9E0); //(CSW_FSCM_OFF & FRC);             //1,8 MIPS
//_FWDT(WDT_ON & WDTPSA_512 & WDTPSB_4);     //WDT = 4096 ms, Prescaler A = 512, B = 4  
_FWDT(WDT_OFF);                                //WDT desligado
_FBORPOR(MCLR_EN & PBOR_OFF & PWRT_OFF);
_FGS(CODE_PROT_OFF); 


//---------------------- Define constantes -------------------------------
// constantes da FFT
#define PONTOS_FFT    128         //numero de pontos da FFT
#define LOG2_ESTAGIOS_BUTTERFLY 7    //numero de estagios "Butterfly" da FFT
#define TAXA_AMOSTRAGEM        7680    //taxa de amostragem do sinal de entrada 
                                        

//----------------- Declara variaveis globais-----------------------------------

fractional sinal_Entrada [128]= {            //vetor que ira armazenar os dados de entrada
                            //valores iniciais declarados para simulacao
                            //e teste do programa
                            //square signal 780 Hz 7.68 kSPS
                            //sinal gerado com o dsPICWorks
    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, 0x8001, 0x8001,
    0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
    0x8001, 0x8001, 0x8001, 0x8001, 0x8001,
    0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF, 0x7FFF,
    0x8001, 0x8001, 0x8001};


unsigned char *TXPtr;
unsigned int n,dado;
float tempFloat, modulo;


int    picoFrequenciaBin = 0;                
unsigned long picoFrequencia = 0;

fractional *temp3_pnt, tempFrac;

// Entrada do sinal para a FFT declarada na memoria Y
fractcomplex sinalComplexo[PONTOS_FFT]__attribute__ ((section (".ydata, data, ymemory"),aligned (PONTOS_FFT * 2 *2))); 

fractional SquareMagnitude [PONTOS_FFT/2];


// ------------- Declara constantes gravadas na flash -------------------------

// declara vetor dos coeficientes Twiddle da FFT na memoria de programa
// Define 128 constantes "twiddle" da FFT: WN (kn) = exp [-(j*2*pi*k*n)/N]
const fractcomplex twiddleFactors[] __attribute__ ((space(prog), aligned (PONTOS_FFT*2)))=
        {
        0x7FFF, 0x0000, 0x7FD9, 0xF9B8, 0x7F62, 0xF374, 0x7E9D, 0xED38,
        0x7D8A, 0xE707, 0x7C2A, 0xE0E6, 0x7A7D, 0xDAD8, 0x7885, 0xD4E1,
        0x7642, 0xCF04, 0x73B6, 0xC946, 0x70E3, 0xC3A9, 0x6DCA, 0xBE32,
        0x6A6E, 0xB8E3, 0x66D0, 0xB3C0, 0x62F2, 0xAECC, 0x5ED7, 0xAA0A,
        0x5A82, 0xA57E, 0x55F6, 0xA129, 0x5134, 0x9D0E, 0x4C40, 0x9930,
        0x471D, 0x9592, 0x41CE, 0x9236, 0x3C57, 0x8F1D, 0x36BA, 0x8C4A,
        0x30FC, 0x89BE, 0x2B1F, 0x877B, 0x2528, 0x8583, 0x1F1A, 0x83D6,
        0x18F9, 0x8276, 0x12C8, 0x8163, 0x0C8C, 0x809E, 0x0648, 0x8027,
        0x0000, 0x8000, 0xF9B8, 0x8027, 0xF374, 0x809E, 0xED38, 0x8163,
        0xE707, 0x8276, 0xE0E6, 0x83D6, 0xDAD8, 0x8583, 0xD4E1, 0x877C,
        0xCF04, 0x89BE, 0xC946, 0x8C4A, 0xC3A9, 0x8F1D, 0xBE32, 0x9236,
        0xB8E3, 0x9592, 0xB3C0, 0x9931, 0xAECC, 0x9D0E, 0xAA0A, 0xA129,
        0xA57E, 0xA57E, 0xA129, 0xAA0A, 0x9D0E, 0xAECC, 0x9931, 0xB3C0,
        0x9592, 0xB8E3, 0x9236, 0xBE32, 0x8F1D, 0xC3A9, 0x8C4A, 0xC946,
        0x89BE, 0xCF04, 0x877C, 0xD4E1, 0x8583, 0xDAD8, 0x83D6, 0xE0E6,
        0x8276, 0xE707, 0x8163, 0xED38, 0x809E, 0xF374, 0x8027, 0xF9B8
        } ;



//----------------- Declara funcoes externas ---------------------------------

// calcula a magnitude quadrada do vetor complexo
extern void SquareMagnitudeCplx(int, fractcomplex*, fractional*); 





//****************************************************************************************
//*********************************** MAIN ***********************************************
//****************************************************************************************


int main (void)            //inicio do programa
{
    int i = 0;
    fractional *p_real = &sinalComplexo[0].real ;
    fractcomplex *p_complexo = &sinalComplexo[0] ;
    fractional *p_fract = &sinalComplexo[0].real ;


    
  while (1)      //loop infinito
  {

     
    p_real = &sinalComplexo[0].real ;              //inicializa ponteiro
    p_complexo = &sinalComplexo[0];                //inicializa ponteiro
    
    /* Move buffer do sinal de entrada (real)para vetor da FFT    */
    for ( i = 0; i < PONTOS_FFT; i++ )             
    {
        *p_real = sinal_Entrada [i] ;        
        *p_real++;                                 //incrementa o ponteiro     
    }

    
    /* Escalona o vetor para o range [-0.5, +0.5]                */
    p_real = &sinalComplexo[0].real ;             //inicializa ponteiro
    for ( i = 0; i < PONTOS_FFT; i++ )         
    {
        *p_real = *p_real >>1 ;                     //desloca 1 bit para a direita
        *p_real++;            
    }                    


    /* Converte vetor real para vetor complexo                    */
    p_real = &sinalComplexo[(PONTOS_FFT/2)-1].real ;    //inicializa o ponteiro para parte real do vetor  
    p_complexo = &sinalComplexo[PONTOS_FFT-1] ; 

    for ( i = PONTOS_FFT; i > 0; i-- )         
    {
        (*p_complexo).real = (*p_real--);    
        (*p_complexo--).imag = 0x0000;            //Não possui parte imaginaria, escreve valor zero
    }


    /* Calcula a FFT - Transformada Rapida de Fourier             */
    FFTComplexIP (LOG2_ESTAGIOS_BUTTERFLY, &sinalComplexo[0], (fractcomplex *) __builtin_psvoffset(&twiddleFactors[0]), (int) __builtin_psvpage(&twiddleFactors[0]));


    /* Organiza dados da saída da FFT - Bit-reverso                */
    BitReverseComplex (LOG2_ESTAGIOS_BUTTERFLY, &sinalComplexo[0]);


    /* Calcula a soma do quadrado das partes real e imaginária 
       transformando vetor complexo em um vetor real 
            Sinal = real^2 + imaginario^2                      */
    SquareMagnitudeCplx(PONTOS_FFT, &sinalComplexo[0], &sinalComplexo[0].real);

    /* Move buffer do vetor complexo para vetor comum             */
    p_fract = &sinalComplexo[0].real ;                //inicia ponteiro
    for ( i = 0; i < PONTOS_FFT/2; i++ )             
    {
        SquareMagnitude [i] = *p_fract;        
        *p_fract++;                                     //incrementa o ponteiro                                               
    }


    /* Encontra a frequency Bin ( indice da frequencia) que possui a maior energia */
    VectorMax(PONTOS_FFT/2, &SquareMagnitude[0], &picoFrequenciaBin);

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

    /* Encontra  o modulo da frequencia de pico                    */
    temp3_pnt = &SquareMagnitude[0] ;      //inicializa ponteiro    
    temp3_pnt = temp3_pnt + picoFrequenciaBin;
    tempFrac = *temp3_pnt;
    tempFloat = tempFrac;
    tempFloat = tempFloat / 32767;            // Converte fracionario para float
    modulo = sqrt (tempFloat);                // Extrai a raiz quadrada da soma dos quadrados da parte real/imaginaria

 
    //tempFloat = Fract2Float (tempFrac);     //converte fracional para float, usando funcao da biblioteca DSP
    //modulo = sqrt (tempFloat);            // Extrai a raiz quadrada da soma dos quadrados da parte real/imaginaria


    // se o programa estiver calculando a FFT corretamente
    // na tabela Watch voce deve encontrar picoFrequencia = 780 Hz
    // de acordo com o sinal que foi gerado para esta simulacao
    
    asm    ("NOP");

   }    
}                        /*  Fim  */




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. dsPIC DSC DSP Algorithm Library. Disponível em: http://www.microchip.com/stellent/idcplg?IdcService=SS_GET_PAGE&nodeId=2680&dDocName=en023598