FFT con un pic?

[hibiscusblau]

Buenas notícias!

Ya he encontrado el fallo, o los fallos! Aquí os los pongo por si alguien le pasa lo mismo que yo:

Correción 1: guardar las muestras en parte real y parte imaginaria

Código:

void __attribute__ ((__interrupt__, no_auto_psv)) _ADCInterrupt(void)
{
    IFS0bits.T3IF = 0;    //Clear the Timer3 Interrupt Flag    
    IFS0bits.ADIF = 0; //Clear the A/D Interrupt flag bit or else the CPU will keep vectoring back to the ISR

	if(tics%2 == 0 ){
		sigCmpx[CountSamples].real = ADCBUF0; 
	}else{
		sigCmpx[CountSamples].imag = ADCBUF0; 
		CountSamples++;
	}
	tics++;
	if (CountSamples >= FFT_BLOCK_LENGTH){	
		doFilterFlag = 1; 
		CountSamples = 0;
		tics = 0;
	}
}

Corrección 2:

Código:
#define FFT_BLOCK_LENGTH	512    /* = Number of frequency points in the FFT */
#define LOG2_BLOCK_LENGTH 	9	/* = Number of "Butterfly" Stages in FFT processing */  2^9 = 512

#define SAMPLING_RATE		8000 	/* = Rate at which input signal was sampled */
                                        /* SAMPLING_RATE is used to calculate the frequency*/
                                        /* of the largest element in the FFT output vector*/

La fft se realiza asi:

char FFTMotor(){
	int i = 0;

	double aux = 0;

	fractional *p_real = &sigCmpx[0].real ; // posa al punter p_real l'@de memòria on es troba sigCmpx[0].real
	fractcomplex *p_cmpx = &sigCmpx[0] ; // posa al punter p_cmpx l'@de memòria on es troba sigCmpx[0]

	IEC0bits.ADIE = 0;		 //paro la interrupció per fer la FFT

	sprintf(uTxt, "\r\n entro  a la FFTMOTOR \r\n" );
	putsUART1((unsigned int *)uTxt);
	ThereIsARockFlag = 0;
	doFilterFlag = 0;

	p_real = &sigCmpx[(FFT_BLOCK_LENGTH/2)-1].real ;	/* Set up pointers to convert real array */  // p_real = &sigCmpx[127].real posa a p_real l'@ de memòria de sigCmpx[127].real	
	p_cmpx = &sigCmpx[FFT_BLOCK_LENGTH-1]; /* to a complex array. The input array initially has all */ // p_cmpx = &sigCmpx[255] posa a p_cmpx l'@ de memòria de sigCmpx[255]
						/* the real input samples followed by a series of zeros */

	for ( i = FFT_BLOCK_LENGTH; i > 0; i-- )
	{	
 		/* Convert the Real input sample array to a Complex input sample array  */			
		(*p_cmpx).real = (*p_real--);
	/* We will simpy zero out the imaginary   part of each data sample */
		(*p_cmpx--).imag = 0x0000;	
	}

	
#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

	sigCmpx[0].real = 0; // remove signal continue

	// 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
	aux = 15.63; //((8000)/FFT_BLOCK_LENGTH);
	peakFrequency = peakFrequencyBin*15.63; //(SAMPLING_RATE/FFT_BLOCK_LENGTH); //aux;

	sprintf(uTxt, "\r\n Frequencia de Pic: %ld \r\n", peakFrequency);
	putsUART1((unsigned int *)uTxt);

	IEC0bits.ADIE = 0;
	return 0;
}

Saludos!

 

[jhefren]

Saludos @dsklast podrias explicar paso a paso como trabaja, en realidad quiero hacer una, pero en assembler, ya que es mi fuerte y no manejo muy bien el lenguaje c. Muchas Gracias

 

[hibiscusblau]

Supongo que quieres que saber el paso a paso de la fft, no?

Si te fijas es un codigo que viene proporcionado por las librerías de microchip, y uso sus funciones.

De hecho, actualmente modifiqué el código y es más simple. Y es el siguiente:

void __attribute__ ((__interrupt__, no_auto_psv)) _ADCInterrupt(void)
{
    //Clear the Timer3 Interrupt Flag
    IFS0bits.T3IF = 0;
    //Clear the A/D Interrupt flag bit or else the CPU will     //keep vectoring back to the ISR
    IFS0bits.ADIF = 0;
sigCmpx[CountSamples].real = ADCBUF0; //>>1;
	sigCmpx[CountSamples].imag = 0; //>>1;
	CountSamples++;

	if( ADCBUF0 > AdData){
		AdData = ADCBUF0;
	} 
	tics++;
	if (CountSamples >= FFT_BLOCK_LENGTH){
		doFFTflag = 1; // CountSamples >= 256, ja ho tinc tot
		CountSamples = 0;
		t++;
		tics = 0;
	}
//	LEDblau = 1;
}

 

[hibiscusblau]

Hola jhefren,

El otro día no pude terminar mi post. El código que te puse es la adquisición de los datos, que si te fijas, guardo los datos en la parte real de sigCmpx, y en la parte imaginaria 0.

sigCmpx[CountSamples].real = ADCBUF0; //>>1;
sigCmpx[CountSamples].imag = 0; //>>1;

Ya que si coges los típicos ejemplos de Microchip, al final mueve las muestras de la parte imaginaria a la mitad para arriba (de FFT_BLOCK_LENGTH/2 a FFT_BLOCK_LENGTH) y en mi caso, la parte de la FFT he sacado esta parte para que el código sea más rendible.


Código de la FFT:

	if(doFFTflag == 1){ //si tengo todas las muestras

	        IEC0bits.ADIE = 0;  //paro la interrupción para hacer la FFT	
		doFFTflag = 0;

		                VectorWindow(FFT_BLOCK_LENGTH,&sigCmpx[0].real,&sigCmpx[0].real,ham_window); 
		fractional *p_real = &sigCmpx[0].real ; // pongo el puntero p_real la @ de memoria donde se encuentra sigCmpx[0].real
		// redimensiono la señal para que no sobrepase el margen que la función FFT necesita para operar correctamente
		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]*/
			*p_real = *p_real >>1 ;		/* So, we shift all data samples by 1 bit to the right. */ 
			*p_real++;			/* Should you desire to optimize this process, perform data scaling when first obtaining the time samples */
					/* Or within the BitReverseComplex function source code */
		}	

                //hago la FFT
		#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

		sigCmpx[0].real = 0; // pongo el primer "BIN" a 0 porque no quiero la tensión continua que tiene a la entrada, (esto no se hasta qué punto es necesario pero por si acaso...)
		
		BitReverseComplex (LOG2_BLOCK_LENGTH, &sigCmpx[0]); // Store output samples in bit-reversed order of their addresses. 
                //gira el array 
		SquareMagnitudeCplx(FFT_BLOCK_LENGTH, &sigCmpx[0], &sigCmpx[0].real); // Compute the square magnitude of the complex FFT output array so we have a Real output vetor
                //calcula la magnitud a real, haciendo la raiz de la parte real^2 y la parte imaginaria^2 de cada vector. O sea, pasa de una función compleja a imaginaria. 
		VectorMax(FFT_BLOCK_LENGTH/2, &sigCmpx[0].real, &peakFrequencyBin);	// Compute the frequency (in Hz) of the largest spectral component 		// Find the frequency Bin ( = index into the sigCmpx[] array) that has the largest energy 	// i.e., the largest spectral component
                  //calcula que "BIN" tiene la potencia más alta.


		peakFrequency = peakFrequencyBin*Freq_for_BIN; 
          //calcula qué frecuencia tiene la potencia más alta. 
       }

Cada BIN es cada casilla del array. O sea, si nuestra frecuencia de muestreo es de 4000muestras/segundo y el número de muestras es 512, la unidad del BIN sera:

BIN = Fs/N = 4000/512 = 7.8125Hz.

O sea, si la frequencia de entrada vale 10Hz, la función Vector_Max, devolvera que peakFrequencyBin es 2, ya que se encuentra en la 2a casilla del array. Cada casilla son 7.8125Hz.

BIN 1 = 1a casilla de 0Hz a 7.8125Hz
BIN 2 = 2a casilla de 7.8125Hz a 15.625Hz
BIn 3 = 3a casilla de 15.625Hz a 23.4375Hz

Por eso, para saber que frequencia corresponde el BIN que devuelve peakFrecuencyBin lo multiplico por Freq_for_BIN que está definido como:

#define Freq_for_BIN ((double) SAMPLINGRATE) / FFT_BLOCK_LENGTH;

Espero que te sirva de ayuda.