The analog I/O functionality is based on the code described in this blog entry :
How To Connect Digital Signal Processing Functions To An ADC And A DAC
We can use the free version of SigLib (no source code) from here : http://www.numerix-dsp.com/free/.
Before proceeding any further, please ensure :
That SigLib is installed and configured correctly for your compiler. Full details are included in the SigLib User's Guide which is provided with the library and also available here : http://www.numerix-dsp.com/docs/.
That GnuPlot and Gnuplot/C are both installed and working correctly : http://realgonegeek.blogspot.co.uk/search?q=gnuplot.
The code uses a double buffering technique so that the analog input code writes the data into a separate array to the one being processed with the FFT functions.
Here is the source code for the top level function (pa_fft.c) :
#include <stdio.h>
#include <math.h>
#include "analog_io.h"
#include <siglib.h> /* SigLib DSP library */
#include <gnuplot_c.h> /* Gnuplot/C */
#if _MSC_VER // If compiler not MSVC then include kbhit.h
#include <conio.h>
#else
#include "kbhit.h"
#endif
#define SAMPLE_RATE 44100
#define FFT_LENGTH ((SLArrayIndex_t)1024)
#define LOG2_FFT_LENGTH ((SLArrayIndex_t)10)
#define HALF_FFT_LENGTH (FFT_LENGTH >> 1)
#define WINDOW_SIZE FFT_LENGTH
/* Declare global variables and arrays */
SLData_t *pRealData, *pImagData, *pWindowCoeffs, *pResults, *pFFTCoeffs;
volatile int Input_SamplesCount;
volatile int Input_Data_Valid_Flag;
volatile int Input_Data_Buffer;
SLData_t *pInputData, *pProcessData;
SLData_t Data1 [FFT_LENGTH], Data2 [FFT_LENGTH];
h_GPC_Plot *h2DPlot; /* Plot object */
void analog_isr (void)
{
// Process channel 0 - Store data for FFT and output zero
*(pInputData+Input_SamplesCount) = (double)adc_in0;
dac_out0 = 0;
// Process channel 1 - just output zero
dac_out1 = 0;
Input_SamplesCount++;
if (Input_SamplesCount == FFT_LENGTH) // If we have filled the buffer, mark data as valid and swap buffers
{
Input_Data_Valid_Flag = 1; // FFT Data is now valid
Input_SamplesCount = 0; // Reset input sampels count
if (Input_Data_Buffer == 1)
{
pInputData = Data2; /* Input to array 1 */
pProcessData = Data1; /* Perform FFT on array 2 */
Input_Data_Buffer = 2;
}
else // (Input_Data_Buffer == 2)
{
pInputData = Data1; /* Input to array 1 */
pProcessData = Data2; /* Perform FFT on array 2 */
Input_Data_Buffer = 1;
}
}
}
void process_FFT (void)
{
SDA_Copy (pProcessData, pRealData, FFT_LENGTH); // Copy data for processing
Input_Data_Valid_Flag = 0; // FFT Data has been copied for processing
/* Apply window to data */
SDA_Window (pRealData, /* Pointer to source array */
pRealData, /* Pointer to destination array */
pWindowCoeffs, /* Pointer to window oefficients */
WINDOW_SIZE); /* Window length */
/* Perform real FFT */
SDA_Rfft (pRealData, /* Pointer to real array */
pImagData, /* Pointer to imaginary array */
pFFTCoeffs, /* Pointer to FFT coefficients */
SIGLIB_NULL_ARRAY_INDEX_PTR, /* Pointer to bit reverse address table - NOT USED */
FFT_LENGTH, /* FFT length */
LOG2_FFT_LENGTH); /* log2 FFT length */
SDA_LogMagnitude (pRealData, /* Pointer to real source array */
pImagData, /* Pointer to imaginary source array */
pResults, /* Pointer to log magnitude destination array */
HALF_FFT_LENGTH); /* Data dataset length */
SDA_Offset (pResults, -138.0, pResults, HALF_FFT_LENGTH); // Offset graph results for 0 dB
gpc_plot_2d (h2DPlot, /* Graph handle */
pResults, /* Dataset */
HALF_FFT_LENGTH, /* Dataset length */
"FFT of ADC Input Data", /* Dataset title */
SIGLIB_ZERO, /* Minimum X value */
(double)(HALF_FFT_LENGTH - 1), /* Maximum X value */
"lines", /* Graph type */
"blue", /* Colour */
GPC_NEW); /* New graph */
}
int main(void)
{
int Error;
Input_SamplesCount = 0; // Initialize flags
Input_Data_Valid_Flag = 0;
Input_Data_Buffer = 0;
pInputData = Data1; /* Input to array 1 */
pProcessData = Data2; /* Perform FFT on array 2 */
/* Allocate memory */
pRealData = SUF_VectorArrayAllocate (FFT_LENGTH);
pImagData = SUF_VectorArrayAllocate (FFT_LENGTH);
pFFTCoeffs = SUF_FftCoefficientAllocate (FFT_LENGTH);
pResults = SUF_VectorArrayAllocate (FFT_LENGTH); /* RMS result array */
pWindowCoeffs = SUF_VectorArrayAllocate (WINDOW_SIZE); /* Window array */
h2DPlot = /* Initialize plot */
gpc_init_2d ("Fast Fourier Transform", /* Plot title */
"Frequency", /* X-Axis label */
"Magnitude", /* Y-Axis label */
120.0, /* Scaling mode */
GPC_NEGATIVE, /* Sign mode */
GPC_MULTIPLOT, /* Multiplot / fast plot mode */
GPC_KEY_DISABLE); /* Legend / key mode */
if (h2DPlot == NULL)
{
printf ("\nPlot creation failure.\n");
exit (1);
}
SIF_Window (pWindowCoeffs, /* Pointer to window oefficient */
SIGLIB_HANNING, /* Window type */
SIGLIB_ZERO, /* Window coefficient */
FFT_LENGTH); /* Window length */
/* Initialise FFT */
SIF_Fft (pFFTCoeffs, /* Pointer to FFT coefficients */
SIGLIB_NULL_ARRAY_INDEX_PTR, /* Pointer to bit reverse address table - NOT USED */
FFT_LENGTH); /* FFT length */
Error = analog_open (SAMPLE_RATE, analog_isr); // Open the analog interface
if (Error == -1)
return 1;
printf("Hit ENTER to stop program.\n");
while (!kbhit()) // Wait until key hit
{
if (Input_Data_Valid_Flag == 1) // If data is valid process the FFT
{
process_FFT ();
}
}
Error = analog_close (); // Close the analog interface
if (Error == -1)
return 1;
gpc_close (h2DPlot);
SUF_MemoryFree (pRealData); /* Free memory */
SUF_MemoryFree (pImagData);
SUF_MemoryFree (pResults);
SUF_MemoryFree (pWindowCoeffs);
SUF_MemoryFree (pFFTCoeffs);
return 0;
}
Notes :
The global variables are decalared as volatile. This prevents the compiler optimizing away the code, which often happens when variables are used in interrupt service routines because the compiler has no way of knowing how the ISR is going to be called.
GCC and some other compilers do not implement the kbhit () function so I have included the code from here : http://cboard.cprogramming.com/linux-programming/51531-faq-cached-input-mygetch.html#post357655.
Testing
Under Microsoft Visual Studio the code can be built using the following command (assuming you have already got a working PortAudio installation) :
cl pa_fft.c analog_io.c -W4 -DPA_USE_ASIO=1 -D "SIGLIB_STATIC_LIB=1" -D "_CRT_SECURE_NO_WARNINGS=1" portaudio_x86.lib siglib.lib gnuplot_c.lib nhl.lib
The XMOS xCORE processor on this board is more than capable of handling the signal processing code shown in this blog but for the purposes of the demonstration I am running the DSP code on the host.
If you run the code yourself you will see the the FFT output plotted in Gnuplot and should look something like this - if you apply a noisy sinewave to the input ;-) :
Stay tuned for a future blog post where I will implement this code and the IIR filter code on an XMOS xCORE :-)
If you have found this solution useful then please do hit the Google (+1) button so that others may be able to find it as well.
Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/
No comments:
Post a Comment