• Welcome to our site! Electro Tech is an online community (with over 170,000 members) who enjoy talking about and building electronic circuits, projects and gadgets. To participate you need to register. Registration is free. Click here to register now.

FIR/FFT/IIR on dsPIC

Status
Not open for further replies.

kchriste

New Member
Forum Supporter
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.
 

krazy

New Member
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.
 

kchriste

New Member
Forum Supporter
Here is some messy code I wrote that works with a GLCD. My original post on that topic is here:
http://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.

Latest threads

EE World Online Articles

Loading
Top