PIC32: How to use DSP Library?

I am using a Microstick II Board with the included PIC32MX250F128B and I would like to perform real-time audio processing. I use an analog input to get the input sound and a PWM output, in combination with a low-pass filter followed by an amplifier, as my audio output (connected to headphones).

I would like to do is to perform more advanced signal processing, like filtering, etc. Microchip website says there is a DSP Library dedicated to that. The lib manual can be found here, but there isn't a lot of information in it, just some function prototypes...

In a first time, I would like to perform low-pass and bandpass filtering. Until now, I managed to perform fft using the following code:

// include header files
#include <plib.h>
#include <p32xxxx.h>
#include <dsplib_dsp.h>
#include <fftc.h>

// Config Bits
#pragma config FNOSC = FRCPLL       // Internal Fast RC oscillator (8 MHz) w/ PLL
#pragma config FPLLIDIV = DIV_2     // Divide FRC before PLL (now 4 MHz)
#pragma config FPLLMUL = MUL_20     // PLL Multiply (now 80 MHz)
#pragma config FPLLODIV = DIV_2     // Divide After PLL (now 40 MHz)
#pragma config FWDTEN = OFF         // Watchdog Timer Disabled
#pragma config ICESEL = ICS_PGx1    // ICE/ICD Comm Channel Select (pins 4,5)
#pragma config JTAGEN = OFF         // Disable JTAG
#pragma config FSOSCEN = OFF        // Disable Secondary Oscillator

// Defines

#define fftc fft16c64 // from fftc.h, for N = 64

#define SYSCLK (40000000L)

#define SAMPLES 64

#define PWM_FREQ   48000            // Output PWM frequency
#define DUTY_CYCLE  1

int tab[SAMPLES];
int i = 0;

{

while(AD1CON1bits.SAMP);        // wait until acquisition is done
while(!AD1CON1bits.DONE);       // wait until conversion done

}

{

AD1CON1 = 0x00E0;       // internal counter ends sampling and starts conversion (manual sample)

// Found on the web (todo: check the datasheet)
}

int main( void)
{
SYSTEMConfigPerformance(SYSCLK);

// Set OC1 to pin 2 with peripheral pin select
RPA0Rbits.RPA0R = 0x0005;

// Configure standard PWM mode for output compare module 1
OC1CON = 0x0006;

for(i = 0; i<SAMPLES; i++)
tab[i] = 0;

// From datasheet:
// PR = [FPB / (PWM Frequency * TMR Prescale Value)] – 1
PR2 = (SYSCLK / PWM_FREQ) - 1;

// Initial duty cycle value
OC1RS = (PR2 + 1) * ((float)DUTY_CYCLE / 100);

T2CONSET = 0x8000;      // Enable Timer2, prescaler 1:1
OC1CONSET = 0x8000;     // Enable Output Compare Module 1

ANSELBbits.ANSB3 = 1;   // set RB3 (AN5) to analog
TRISBbits.TRISB3 = 1;   // set RB3 as an input
TRISBbits.TRISB5 = 0;   // set RB5 as an output (note RB5 is a digital only pin)

int pos=0, dat = 0;

int log2N = 6; // log2(64) = 6
int N = 1 << log2N; // N = 2^6 = 64
int din[N];
int dout[N];
int dout2[N];
int scratch[N];

while(1)
{
//foo = analogRead 5); // note that we call pin AN5 (RB3) by it's analog number

//dat += din[pos]*3/10;
//if(dat > 1023) dat -= 512;

din[pos] =  dat;

mips_fft16(dout, din, fftc, scratch, 1);
mips_fft16(dout2, dout, fftc, scratch, 1);

if(++pos >= SAMPLES) pos = 0;

OC1RS = (PR2 + 1) * ( ((float)dout2[pos])/1023); // Write new duty cycle
}

return 0;
}


First of all, I would like to know if what I am doing is correct. I mean, there is no ifft function, does the second fft perform an ifft ? Secondly, by doing this, I can hear a high frequency sound, which may result from latency that delays the PWM duty cycle change OC1RS. How can I fix this ? I saw people usually use Timer2 to change the duty cycle, but each time I tried, I get nothing at the output (no sound). How can I implement this in my case ?

• Do you have a copy of "Numerical Recipes?" It is a very good book on these and other subjects. Including doing filtering with FFTs/IFFTs. But it's not hard to see how to do it with brute force. Make sure you understand windowing functions and their purpose, too. E Brigham's "Fast Fourier Transform and Its Applications" in pretty much ANY edition is a top-notch book on the subject, generally. And your question seems pretty general, so that's also a good book to get and read thoroughly. – jonk Mar 22 '14 at 19:44
• I have heard about it, thanks for this good suggestion. But my question is more about "how to do it correctly", I mean in term of performance. I know how to use FFT, I used them a lot with Matlab, but never with a PIC or dsPIC. So, I would like to use DSP lib function as relevant as possible. – Spin Mar 22 '14 at 20:06
• You need chapter 3 of ww1.microchip.com/downloads/en/DeviceDoc/51685E.pdf – Leon Heller Mar 22 '14 at 20:38
• Performance is at least four questions: What can "canned" libraries achieve? How can I change the problem statement to make better use of "canned" libraries for performance? What could I achieve if I wrote custom code for my needs? Do any or all of these meet my requirements? Since you know your application better than any of us, and since the above are really "application questions," I'm not sure how we help much here without a lot more information about the application details. (Low pass is easy, though. What have you tried?) – jonk Mar 22 '14 at 20:43
• Thanks to your comments, I was able to perform FFT. But I still have some questions about the results I get, and particularly about the use of Timer2, which could eventually fix the latency problem. – Spin Mar 23 '14 at 16:02

I find this best learned using the DCMI tool and simulation, although mipsfir doesn't sim for some reason.

mips_fft16/32 does fft and ifft. Its all the same, really. Transfers input from one domain to the other.

Your implementation has flaws. The output of the mipsfft is a complex number, I am surprised that even compiled. Do a study on complex numbers. I don't see your fftc define, but a size of 1 doesn't make sense in the mips_fft16 call. Its really pretty simple. If you want to create a sine wave in a sample buffer of 32, you put the (amplitute*32)value in din[1].re. din[2].re for a double sine, and so on. So for a sample rate of 16000, 16000/32 = 500 * 1 = 500hz. Note that your output is in dout[x].re. So a call to fft/ifft will look more like this for a 32 sized fft. mips_fft16(dout, din, fftc16c32, scratch, 6);

Also, remember that the output of an fft is going to be real and imaginary, so if you just want peaks don't forget that A squared + B squared = C squared.

Real and imaginary = phase of output, greatly simplified.

OC1RS = (PR2 + 1) * ( ((float)dout2[pos])/1023); // Write new duty cycle


This isn't the proper way to do this at all. The OC1RS update can be anywhere, but I suggest putting in inside a timer at some chosen sample rate, say 16000. Then feed OC1RS with an incrementing data pointer of some sort to walk through your output.