FIR/FFT/IIR on dsPIC

Status
Not open for further replies.

krazy

New Member
Does *anyone* here have any experience with performing any kind of filtering on the dsPIC?
 
Microchip has a DSP library that does all those functions. I've used their FFT routines with success. Haven't tried their FIR and IIR routines.
 
Do you have any sample code around that I could take a look at? I cannot get any filtering of anything going for the life of me.
 
Here is some messy code I wrote that works with a GLCD. My original post on that topic is here:
https://www.electro-tech-online.com/threads/dspic-glcd-spectrum-analyzer-d.85530/

Code:
#include <p33FJ64MC802.h>
#include <dsp.h>
#include "fft.h"
#include "sqrtX.h"


_FOSC(OSCIOFNC_ON & FCKSM_CSDCMD & POSCMD_NONE);	//Oscillator Configuration
_FOSCSEL(FNOSC_FRC);								//Oscillator Selection
_FWDT(FWDTEN_OFF);									//Turn off WatchDog Timer
_FGS(GCP_OFF);										//Turn off code protect
_FPOR( FPWRT_PWR1 );								//Turn off power up timer


#define CS2		LATBbits.LATB12
#define CS1		LATBbits.LATB11
#define	DI		LATBbits.LATB10
#define RW		LATBbits.LATB9
#define E		LATBbits.LATB8
#define LCD_RESETTING	PORTBbits.RB4

void e_togg(void)
{
	unsigned int x;
	E = 1;
	for( x=0; x<6; x++){;} //Delay apx 10us
	E = 0;
}

void cmd_write(unsigned char data)
{
	unsigned int temp;
	DI = 0;
	RW = 0;
	temp = LATB;
	temp &= 0xff00;
	temp |= data;
	LATB = temp;
	e_togg();
}

void clearGLCD(void)
{
	unsigned int page, y;
	for( page=0xB8; page<0xC0; page++)
	{
		cmd_write(0x40);		//Set Y address to 0
		cmd_write(page);		//Set pages 0-7
			   //FEDCBA9876543210
		LATB = 0b0001110000000000u;	//CS1,CS2 & DI=1	All others & data =0
		for( y=0; y<64; y++)
			{ e_togg(); }		//Clear byte
	}	
}

void WaitTilReady(void)
{
	unsigned int x;
	TRISB = 0x00ff;
	LATB = 0b0000101000000000u; // CS1 = 1; CS2 = 0; RW = 1; DI = 0;
	for( x=0; x<72; x++){;} //Delay apx 100us
	E = 1;
	for( x=0; x<72; x++){;} //Delay apx 100us
	while((PORTB & 0x0090) != 0);
	RW = 0;
	E = 0;
	for( x=0; x<72; x++){;} //Delay apx 100us
	CS1 = 0;
	CS2 = 1;
	RW = 1;
	DI = 0;
	for( x=0; x<72; x++){;} //Delay apx 100us
	E = 1;
	for( x=0; x<72; x++){;} //Delay apx 100us
	while((PORTB & 0x0090) != 0);
//	while(LCD_RESETTING);
	RW = 0;
	E = 0;
	TRISB = 0x0000;
}



unsigned char sqrtI( unsigned int sqrtArg )
	{
	unsigned char answer, x;
	unsigned int temp;
	if ( sqrtArg < 0x2 ) { return sqrtArg; }
	answer = 0x00;
	for( x=0x80; x>0; x>>=1 )
		{
		answer |= x;
		temp = answer * answer;
		if (temp > sqrtArg)
			{ answer ^= x; }
		if (temp == sqrtArg)
			{ break; }
		}
	return answer;
	}








/* Global Definitions */
unsigned int x,page,RevPage,temp,bin;
unsigned int *TestPtr;
fractcomplex sigCmpx[FFT_BLOCK_LENGTH] __attribute__ ((space(ymemory),far,aligned(FFT_BLOCK_LENGTH * 2 *2))); 

#ifndef FFTTWIDCOEFFS_IN_PROGMEM
fractcomplex twiddleFactors[FFT_BLOCK_LENGTH/2] 	/* Declare Twiddle Factor array in X-space*/
__attribute__ ((section (".xbss, bss, xmemory"), aligned (FFT_BLOCK_LENGTH*2)));
#else
extern const fractcomplex twiddleFactors[FFT_BLOCK_LENGTH/2]	/* Twiddle Factor array in Program memory */
__attribute__ ((space(auto_psv), aligned (FFT_BLOCK_LENGTH*2)));
#endif

void initAdc1(void);

int	peakFrequencyBin = 0;				/* Declare post-FFT variables to compute the */
unsigned long peakFrequency = 0;			/* frequency of the largest spectral component */

int main(void)
{
	int i = 0;
	fractional *p_real = &sigCmpx[0].real ;
	fractcomplex *p_cmpx = &sigCmpx[0] ;

	AD1PCFGL = 0xfffe;	//make ADC pins all digital except RA0
	TRISA = 0x0001;			//Make all PORTs all outputs except AN0
	TRISB = 0;
	initAdc1();
	WaitTilReady();
	LATB = 0b0001110000000000u;	//CS1,CS2 & DI=1	All others & data =0
	for ( x=0; x<3500; x++){;}
	cmd_write(0x3f);		//turn on display
	clearGLCD();

	for ( i=0; i<256; i++)
		{ x = sqrtX(i);}
	x=0;
	for ( i=0; i<256; i++)
		{ x = sqrtI(i);}
	x=0;
while(1){
//	for ( x=0; x<35000; x++){;}
#ifndef FFTTWIDCOEFFS_IN_PROGMEM					/* Generate TwiddleFactor Coefficients */
	TwidFactorInit (LOG2_BLOCK_LENGTH, &twiddleFactors[0], 0);	/* We need to do this only once at start-up */
#endif
	for ( i = 0; i < FFT_BLOCK_LENGTH; i++ )/* The FFT function requires input data */
	{					/* to be in the fractional fixed-point range [-0.5, +0.5]*/

	LATBbits.LATB15 = 1; //ADC Sample flag 1.7uS
		AD1CON1bits.DONE = 0;
		AD1CON1bits.SAMP = 1; 		// start sampling then... after 31Tad go to conversion.
		while (!AD1CON1bits.DONE);	// conversion done?
	LATBbits.LATB15 = 0; //ADC Sample flag 200ms
		sigCmpx[i].real = (signed)ADC1BUF0 >> 1;	/* So, we shift all data samples by 1 bit to the right. */
		sigCmpx[i].imag = 0x0000;
	}

	/* Perform FFT operation */
#ifndef FFTTWIDCOEFFS_IN_PROGMEM
	FFTComplexIP (LOG2_BLOCK_LENGTH, &sigCmpx[0], &twiddleFactors[0], COEFFS_IN_DATA);
#else
	FFTComplexIP (LOG2_BLOCK_LENGTH, &sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&twiddleFactors[0]), (int) __builtin_psvpage(&twiddleFactors[0]));
#endif

	/* Store output samples in bit-reversed order of their addresses */
	BitReverseComplex (LOG2_BLOCK_LENGTH, &sigCmpx[0]);

	/* Compute the square magnitude of the complex FFT output array so we have a Real output vetor */
	SquareMagnitudeCplx(FFT_BLOCK_LENGTH, &sigCmpx[0], &sigCmpx[0].real);

	/* Find the frequency Bin ( = index into the SigCmpx[] array) that has the largest energy*/
	/* i.e., the largest spectral component */
	VectorMax(FFT_BLOCK_LENGTH/2, &sigCmpx[0].real, &peakFrequencyBin);

	/* Compute the frequency (in Hz) of the largest spectral component */
	peakFrequency = peakFrequencyBin*(SAMPLING_RATE/FFT_BLOCK_LENGTH);



	for( x=0; x<465; x++){;} //Delay apx 1ms
	E = 0;
	CS1 = 1;
	CS2 = 1;
	for( x=0; x<72; x++){;} //Delay apx 100us
	cmd_write(0x40);		//Set Y address to 0
	cmd_write(0xb8);		//go to page 0

	for( x=0; x<465; x++){;} //Delay apx 1ms
	DI = 1;
	RW = 0;

	LATBbits.LATB13 = 1; //Display write flag 150ms

	TestPtr = (unsigned int *)&sigCmpx[0];
	RevPage = 8;
	for( page = 0xBF; page >= 0xB8; page-- )
		{
		CS1 = 1; CS2 = 0; cmd_write(page); cmd_write(0x40);
		for( bin = 0; bin <= 127; bin++ )
			{
			if( bin == 64)
				{ CS1 = 0; CS2 = 1; cmd_write(page); cmd_write(0x40); }
			DI = 1; RW = 0;
			if ( TestPtr[bin] > 0x0FFF ) { x = 0x0FFF; }
			else { x = TestPtr[bin] << 3; }
 			if ( x > 0x0FFF ) { x = 0x0FFF; }
			temp = sqrtI(x); //Calculate linear amplitude of frequency bin.
//			temp = sqrtI(TestPtr[bin]); //Calculate linear amplitude of frequency bin.
			if( temp >= RevPage )
				{ LATB |= 0x00FF; } // Make 8 bit bargraph fully ON
			else
				{
				temp = 255 << (RevPage-temp); // Convert to 8 bit bargraph
				temp &= 0x00ff;
				LATB &= 0xff00;
				LATB |= temp;
				}
			e_togg();
			}
		RevPage += 8;
		}	 
	LATBbits.LATB13 = 0; //Display write flag 150ms


}

	return 0;
}









/*=============================================================================
initAdc1() is used to configure A/D to convert channel 5 on Timer event. 
It generates event to DMA on every sample/convert sequence. ADC clock is configured at 625Khz 
=============================================================================*/
void initAdc1(void)
{

		AD1CON1 = 0x00E0; 			// SSRC bit = 111 implies internal counter ends sampling and starts converting.
		AD1CON1bits.FORM   = 3;		// Data Output Format: Signed Fraction (Q15 format)
		AD1CON1bits.AD12B  = 1;		// 12-bit ADC operation

		AD1CON2 = 0;				// Converts CH0
		AD1CON3bits.SAMC = 4;	 	// Sample time = 4Tad, Tad = internal 2Tcy
 		AD1CON3bits.ADRC=0;			// ADC Clock is derived from Systems Clock
		AD1CON3bits.ADCS = 8;		// ADC Conversion Clock Tad=Tcy*(ADCS+1)= (1/7.37M)*64 = 8.7us (115Khz)
									// ADC Conversion Time for 12-bit Tc=14*Tad = 22.4us 

        AD1CHS0bits.CH0SA=0;		// MUXA +ve input selection (AN0) for CH0
		AD1CHS0bits.CH0NA=0;		// MUXA -ve input selection (Vref-) for CH0

  		AD1PCFGL=0xFFFE;			// AN0 as Analog Input
        
        IFS0bits.AD1IF = 0;			// Clear the A/D interrupt flag bit
        IEC0bits.AD1IE = 0;			// Do Not Enable A/D interrupt 
        AD1CON1bits.ADON = 1;		// Turn on the A/D converter	

       
}
Like I said the code is messy and contains redundant code such as two versions of square root for testing purposes. The twiddleFactors table comes with the MicroChip DSP library, so I won't include it here.

Code for the assembler version of square root:
Code:
;-----------------------------------------------------------------------;
; sqrtX.s: Integer square root function.
;
; Input:
;
;       (w0) unsigned int x
;
; Output:
;
;       (w0) unsigned char square root of x
;
; Description:
;
;       Computes the square root of the argument x.
;       
;-----------------------------------------------------------------------;

        .global _sqrtX

_sqrtX:
 		cp		w0,#2
		bra		geu, _domore
		return
_domore:
		clr		w1				;answer reg
		mov		#0x0080,w2		;loop reg x
_sqrtLoop:
		ior		w2,w1,w1		;answer |= x;
		mul.uu	w1,w1,w4		;temp w4 = answer * answer;
		cp		w0,w4			;if (temp w4 > sqrtArg)
		bra		geu, _donext
		xor		w2,w1,w1		;answer ^= x;
_donext:
		cp		w0,w4
		bra		z, _quit		;if (temp w4 == sqrtArg)
		lsr		w2,w2
		bra		nz,_sqrtLoop
_quit:
		mov		w1,w0
        return                  ; Done
;-----------------------------------------------------------------------;
        .end
 
Status
Not open for further replies.
Cookies are required to use this site. You must accept them to continue using the site. Learn more…