Tuesday, August 21, 2012

Remember remember, the... whenever your last DSP class was

After ensuring that I had a good, valid stream of data coming in, I figured it was time to brush off the DSP textbook that had remained undisturbed on my shelf since college and start doing some analysis on the input audio signal.  I started by reading through the documentation and fft example that CMSIS ships with.  The fft example showed how to calculate a 1024 point transform on a complex input signal, which was not exactly what I wanted to do, but was close.  I added the example code to my project and verified that I was able to properly link to the CMSIS .lib that I had previously created when I wrote my app note on the topic (found at http://www.ti.com/lit/an/spma041a/spma041a.pdf).  That done, I switched over to cygwin and wrote a quick c program to emulate what the ADC would see if it was observing a perfect 440 Hz sine wave:


#include <stdio.h>
#include <math.h>

#define PI 3.141592653589793
#define ARR_SIZE 2048
#define SAMPLE_FREQ 44600
#define FREQ 440

int main(int argc, char** argv)
{
        float period = 1.0/SAMPLE_FREQ;
        float time = 0;
        float freqRadians = 2*PI*FREQ;
        unsigned long x[ARR_SIZE*2];
        int i, j = 0;
        printf("const float ti_sample_sine_vector[%d] = {\n", ARR_SIZE);
        for(i=0; i< ARR_SIZE; i++)
        {
                j++;
                //
                // ADC sees signal centered at 0x800 w/ 0x280 max swing
                //
                x[i] = (unsigned long)(640 * sin(freqRadians*time) + 0x800);
                printf("0x%04x,\t", x[i]);
                if(j == 8)
                {
                        printf("\n");
                        j = 0;
                }
                time += period;
        }

        printf("};\n\n\n");
        return 0;
}

Next, I wrote a similar snippet of code to generate the parameters of a windowing function.  I had read in a number of places that an fft is designed to work on a continuous, infinite signal, so when you instead use it for frequency analysis on a finite, non-continuous signal, it is best practice to multiply the samples by a window to reduce the amount of noise across frequency bins.  I think.  I'm nowhere close to a DSP expert, but this sounded like good information to me, so I wrote another snippet to generate a vector that would represent a hamming window:

#include <stdio.h>
#include <math.h>

#define PI 3.141592653589793
#define ARR_SIZE 2048

int main(int argc, char** argv)
{

        double x[ARR_SIZE*2];
        int i, j = 0;
        printf("const float ti_hamming_window_vector[%d] = {\n", ARR_SIZE);
        for(i=0; i< ARR_SIZE; i++)
        {
                j++;
                x[i] = 0.54 - 0.46*cos(2*PI*i/ARR_SIZE);
                //printf("%.18f,\t%.18f,\t", x[i], x[i]);
                printf("%.18f,\t", x[i]);
                //if(j == 4)
                if(j == 8)
                {
                        printf("\n");
                        j = 0;
                }
        }

        printf("};\n\n\n");
        return 0;
}

I copied the output of these two programs into two C files, and used them to write static vectors into the lm4f's flash.  Once that was done, I wrote a bit of code in my launchpad project that would multiply the sample by the hamming window, take a 2048 point real fft of that sample (which is implemented in CMSIS by doing a 1024 point complex fft, then running the result through a split fft process.  I don't know exactly what that means, but I think the result is a real fft of 2048 points).  I then took the complex magnitude of the results, found the maximum value in the result array, determined what frequency range that bin corresponded to, and used my debug UART to print that frequency range.  I had to massage the data and play around a bit with the input flags, but eventually I saw what I had been striving for: a message stating that the maximum frequency component of the signal I ran my DSP function on was between 435 and 458 Hz!

#define NUM_SAMPLES 2048
#define SAMPLING_FREQ 44600
#define INVERT_FFT 0
#define BIT_ORDER_FFT 1

extern q15_t ti_sample_sine_vector[NUM_SAMPLES];
float32_t nextTest[NUM_SAMPLES];

unsigned long g_ulADCValues[NUM_SAMPLES];

float32_t g_fFFTResult[NUM_SAMPLES * 2];


arm_rfft_instance_f32 fftStructure;
arm_cfft_radix4_instance_f32 cfftStructure;

float HzPerBin;

void
SampleSineTest(void)
{
unsigned long i;

arm_rfft_init_f32(&fftStructure, &cfftStructure, NUM_SAMPLES, INVERT_FFT, BIT_ORDER_FFT);
HzPerBin = (float)SAMPLING_FREQ / (float)NUM_SAMPLES;

for(i=0;i<NUM_SAMPLES;i++)
{
//
                // Recenter the signal at 0... might not be necessary
                //
g_fFFTResult[i] = ((float)ti_sample_sine_vector[i] - (float)0x800);
}

UARTprintf("test processing using sample 440 Hz sine wave\n");

arm_mult_f32(g_fFFTResult, ti_hamming_window_vector, g_fFFTResult, NUM_SAMPLES);
UARTprintf("Done with Hamming window multiplication\n");

UARTprintf("Calculating FFT\n");
arm_rfft_f32(&fftStructure, g_fFFTResult, g_fFFTResult);

UARTprintf("Calculating complex magnitude of each frequency bin\n");
arm_cmplx_mag_f32(g_fFFTResult, g_fFFTResult, NUM_SAMPLES * 2);


UARTprintf("Finding largest frequency bin\n");

arm_max_f32(g_fFFTResult, NUM_SAMPLES, &maxValue, &i);


UARTprintf("Peak is between %06d and %06d Hz: %06d\n\n", (int)(HzPerBin*i), (int)(HzPerBin*(i+1)), (int)maxValue);
}


My next step was to very, very anxiously put everything together...

No comments:

Post a Comment