Sunday, 30 August 2015

Why Use A High Level Language For DSP ?

The field of Digital Signal Processing is constantly pushing the price / performance envelope of technology and traditionally this has required systems developers to use assembly language for the majority of the time critical signal processing routines. Today's commercial pressures have moved the "goal-posts" dramatically and typical project development timescales require a larger part of the application to be developed using a high level language. Another benefit to using a high level language for the system development is that a system can be rapidly prototyped to prove the algorithms and then hand optimised using assembly code for the time critical areas.


Primary Reasons For Using High Level Languages

  • High productivity
  • Portability
  • Maintainability
  • Code reuse
  • Optimising system cost / performance
  • Rapid prototyping and algorithm proving
  • Integration with real-time kernels and operating systems
  • Ease of debug
  • Availability of algorithms


The latest generation of compilers allows high level code to be compiled to a quality of assembly code that is very close to that which would be generated by hand. The development process is therefore very much easier than writing the algorithm in assembly code from scratch. An increasingly common development route is to develop the algorithms on a PC or Workstation and then rewrite the application for the target processor. Using the same language for development and deployment often allows the same code to be used for both, with the different I/O requirements handled through the use of conditional compilation of the source.


Modern high performance DSPs are also changing the way we view algorithmic efficiency and an increasing number of projects are written in a high level language because the savings at development time are far greater than the extra cost overhead of using faster processors at deployment.The architectures of the latest DSPs are also becoming more complex, for example with the integration of parallel execution units. This means that it is increasingly difficult for programmers to learn how to fully optimise their algorithms. When the complexity issue is coupled with the fact that the majority of DSP algorithms are block oriented vector processing algorithms and it is now becoming possible for high level language compilers to produce code that is 100% optimised.


DSP Tech Brief : The Zoom-FFT

The Zoom-FFT is a process where an input signal is mixed down to baseband and then decimated, prior to passing it into a standard FFT. The advantage is for example that if you have a sample rate of 10 MHz and require at least 10Hz resolution over a small frequency band (say 1 KHz) then you do not need a 1 Mega point FFT, just decimate by a factor of 4096 and use a 256 point FFT which is obviously quicker.

 Advantages of the Zoom FFT are :

  • Increased frequency domain resolution
  • Reduced hardware cost and complexity
  • Wider spectral range

Applications of the Zoom FFT include :

  • Ultrasonic blood flow analysis
  • R.F. communications
  • Mechanical stress analysis
  • Doppler radar

The following diagram shows the zoom process :



While the following diagram shows the basic architecture of the Zoom-FFT :



One common question is : Is the zoom FFT the same as the chirp z-transform.

The answer is : Absolutely not. The FFT calculates the FFT at N equally spaced points around the unit circle in the z-plane, the chirp z-transform modifies the locations of these points along a contour that can lie anywhere on the z-plane. In contrast, the zoom-FFT uses digital down conversion techniques to localise the standard FFT to a narrow band of frequencies that are centered on a higher frequency. The chirp z-transform is often used to analyze signals such as speech, that have certain frequency domain charactgeristics. The zoom-FFT is used to reduce the sample rate required when analysing narrowband signals - E.G. in HF communications.

These functions, and more, are available in the SigLib DSP Library.

DSP Tech Brief : Signal Scrambling And Linear Feedback Shift Registers (LFSRs)

Modulated signals are often scrambled when they are modulated. There are many reasons for scrambling the signal, including : spreading the spectrum and reducing the correlation between separate channels. Most scramblers utilise a Maximum Length Pseudo-Random Binary Sequence (MLPRBS or MLS), this uses a shift register an exclusive-or gates to randomise the data, as shown in the following diagrams. There are two general arrangements for implementing a scrambler, the non-synchronising and the self-synchronising types, each has its own benefits and draw backs, the non-synchronising type requires that the transmitter and receiver are synchronised before data transmission starts however once synchronised an input bit error to the receive scrambler will cause a single bit error in the output. The self-synchronising scrambler(shown below), as the name suggests will automatically synchronise to the received bit stream but a single receive bit error will cause N+1 bit errors in the output stream, where N is the length of the shift register.



This algorithm is included in the SigLib DSP Library.

DSP Tech Brief : One-Pole Filters

One-pole filters are very useful and efficient techniques for filtering signals. The feedback a value can also be used to design efficient peak-hold filters in the time-domain and are also used in power spectrum analysers and can help extract features of the signals that we are analyzing.



The following equation defines the one-pole filter :

y(n) = x(n) - a * y(n-1), 0 < a < 1.0

The single coefficient defines the gain and time-constant of the filter as follows :

Time constant :



Gain :



This algorithm is included in the SigLib DSP Library.

DSP Tech Brief : Adaptive Filters

The most common adaptive filters are of the FIR type, particularly the Least Mean Square (LMS) adaptive filter.

The LMS filter configuration and adaptation equations are shown in the following diagrams :



The following diagram shows a plot of X^2 gainst Y^2 i.e. the adaptation surface of a 2 tap filter :



This algorithm is included in the SigLib DSP Library.

DSP Tech Brief : The Goertzel Algorithm

Many applications require the detection of a few discrete sinusoids. The Goertzel filter is an IIR filter that uses the feedback to generate a very high Q bandpass filter where the coefficients are easily generated from the required centre frequency, according to the following equations. The most common configuration for using this technique is to measure the signal energy before and after the filter and to compare the two. If the energies are similar then the input signal is centred in the pass-band, if the output energy is significantly lower than the input energy then the signal is outside the pass band. The Goertzel algorithm is most commonly implemented as a second order recursive IIR filter, as shown below.



This algorithm is included in the SigLib DSP Library.

DSP Tech Brief : Envelope Detection

A common and very efficient technique for envelope detection is based on the Hilbert Transform, (a 90 degree phase shifter) and summing the phase shifted and original signals. The Hilbert transform is typically implemented as an FIR filter so the original signal must be delayed to match the group delay of the Hilbert transform. This process can be followed by absolute and then peak hold functions, the latter is often implemented as a one pole IIR filter.



The following diagrams give an idea of how this envelope approximation process works.



The following equation shows how the coefficients of the FIR filter implementation of the Hilbert transform are calculated.



The one pole IIR filter is defined by the following equation :

y(n) = x(n) + a * y(n-1), where 0 < a < 1.0

Another option for detecting the envelope is to use the square root of the energies of the original and the Hilbert transformed signals, as follows :

Envelope = sqrt(Hilbert^2 + signal^2)

In general, this will give similar results to the absolute value technique but can be more run-time efficient.

This algorithm is included in the SigLib DSP Library.


The Next Round Of The University Of Oxford, UK DSP Courses Take Place In January 2016

As part of the University Of Oxford Summer Engineering Program for Industry, the DSP courses are returning in July. Here is a summary of the two courses.

Digital Signal Processing (theory and application) - Monday 25th to Wednesday 27th January 2016
http://www.conted.ox.ac.uk/courses/details.php?id=H600-24
This course provides a good understanding of DSP principles and their implementation and equips the delegate to put the ideas into practice and/or to tackle more advanced aspects of DSP. 'Hands-on' laboratory sessions are interspersed with the lectures to illustrate the taught material and allow you to pursue your own areas of interest in DSP. The hands-on sessions use specially written software running on PCs.

Subjects include :

  • Theoretical Foundations
  • Digital Filtering
  • Fourier Transforms And Frequency Domain Processing
  • DSP Hardware And Programming
  • ASIC Implementation
  • Typical DSP Applications

The course is presented by experts from industry for Engineers in industry.

Digital Signal Processing Implementation (algorithms to optimization) - Thursday 28th January 2016
A one-day supplement to the Digital Signal Processing course that takes the theory and translates it into practice.
http://www.conted.ox.ac.uk/courses/details.php?id=H600-25
The course will include a mixed lecture and demonstration format and has been written to be independent of target processor architecture.
The course will show how to take common DSP algorithms and map them onto common processor architectures. It will also give a guide line for how to choose a DSP device, in particular how to choose and use the correct data word length for any application.

Attendee Feedback From Previous Courses :

  • It was informative, enjoyable and stimulating
  • Excellent content, very lively thanks to the 2 excellent presenters - Anonymous
  • A very good introduction to DSP theory
  • Excellent lecturers! Really useful information and very understandable
  • Great mix of theory and practice
  • The lecturers gave a detailed and excellent explanation of the fundamental topics of DSP with real world engineering practice.
  • This session closes the gap and clears up much confusion between classroom DSP theories and actual DSP implementation.
  • Very good session, with in-depth discussion on the math and background.


The courses will be held at the University of Oxford, UK

Friday, 26 June 2015

Socket API C Programs That Compile For Both UNIX/Linux and Windows

This web page is a good starting point for anyone who wants to learn how to program with the UNIX socket API :
http://www.cs.dartmouth.edu/~campbell/cs60/socketprogramming.html

Many of the IP connected applications I have written need to support both Windows and Unix/Linux. There are some minor differences between the two socket APIs but these are primarily related to configuration.

Once the code has been configured the vast majority of the code will be the same for both operating systems.

These programs use the following predefined macros to select between Windows and Unix/Linux :

#if (_WIN32) || (_WIN64)

Here are versions of the code included in the above URL, that compile on both Windows and Linux.

echoServer.c


#if (_WIN32) || (_WIN64)
#undef UNICODE
#define WIN32_LEAN_AND_MEAN

#include <windows.h>
#include <winsock2.h>
#include <ws2tcpip.h>
#pragma comment (lib, "Ws2_32.lib")
#else
#include <sys/socket.h>
#include <netinet/in.h>
#include <sys/types.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#define MAXLINE 4096 /*max text line length*/
#define SERV_PORT 3000 /*port*/
#define LISTENQ 8 /*maximum number of client connections */

int main (int argc, char **argv)
{
 int n;
 socklen_t clilen;
 char buf[MAXLINE];
#if (_WIN32) || (_WIN64)
 SOCKET listenfd, connfd;
 WSADATA wsaData;
 int iResult;
#else
 int listenfd, connfd;
#endif

 struct sockaddr_in cliaddr, servaddr;

#if (_WIN32) || (_WIN64)
    // Initialize Winsock
    iResult = WSAStartup(MAKEWORD(2,2), &wsaData);
    if (iResult != 0) {
        printf("WSAStartup failed with error: %d\n", iResult);
        return 1;
    }
#endif

 //creation of the socket
 listenfd = socket (AF_INET, SOCK_STREAM, 0);

 //preparation of the socket address 
 servaddr.sin_family = AF_INET;
 servaddr.sin_addr.s_addr = htonl(INADDR_ANY);
 servaddr.sin_port = htons(SERV_PORT);

 bind (listenfd, (struct sockaddr *) &servaddr, sizeof(servaddr));

 listen (listenfd, LISTENQ);

 printf("%s\n","Server running...waiting for connections.");

 for ( ; ; ) {

  clilen = sizeof(cliaddr);
  connfd = accept (listenfd, (struct sockaddr *) &cliaddr, &clilen);
  printf("%s\n","Received request...");

  while ( (n = recv(connfd, buf, MAXLINE,0)) > 0)  {
   printf("%s","String received from and resent to the client:");
   puts(buf);
   send(connfd, buf, n, 0);
  }

 if (n < 0) {
  perror("Read error"); 
  exit(1);
 }
 close(connfd);

 }
 //close listening socket
 close (listenfd);

#if (_WIN32) || (_WIN64)
 WSACleanup();
#endif
}

echoClient.c


#if (_WIN32) || (_WIN64)
#undef UNICODE
#define WIN32_LEAN_AND_MEAN

#include <windows.h>
#include <winsock2.h>
#include <ws2tcpip.h>
#pragma comment (lib, "Ws2_32.lib")
#else
#include <sys/socket.h>
#include <netinet/in.h>
#include <sys/types.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>


#define MAXLINE 4096 /*max text line length*/
#define SERV_PORT 3000 /*port*/

int
main(int argc, char **argv) 
{
 struct sockaddr_in servaddr;
 char sendline[MAXLINE], recvline[MAXLINE];
#if (_WIN32) || (_WIN64)
 SOCKET sockfd;
 WSADATA wsaData;
 int iResult;
#else
 int sockfd;
#endif

#if (_WIN32) || (_WIN64)
    // Initialize Winsock
    iResult = WSAStartup(MAKEWORD(2,2), &wsaData);
    if (iResult != 0) {
        printf("WSAStartup failed with error: %d\n", iResult);
        return 1;
    }
#endif

 //basic check of the arguments
 //additional checks can be inserted
 if (argc !=2) {
  perror("Usage: TCPClient <IP address of the server>"); 
  exit(1);
 }
 //Create a socket for the client
 //If sockfd<0 there was an error in the creation of the socket
 if ((sockfd = socket (AF_INET, SOCK_STREAM, 0)) <0) {
  perror("Problem in creating the socket");
  exit(2);
 }
 //Creation of the socket
 memset(&servaddr, 0, sizeof(servaddr));
 servaddr.sin_family = AF_INET;
 servaddr.sin_addr.s_addr= inet_addr(argv[1]);
 servaddr.sin_port =  htons(SERV_PORT); //convert to big-endian order
 //Connection of the client to the socket 
 if (connect(sockfd, (struct sockaddr *) &servaddr, sizeof(servaddr))<0) {
  perror("Problem in connecting to the server");
  exit(3);
 }
 while (fgets(sendline, MAXLINE, stdin) != NULL) {
  send(sockfd, sendline, strlen(sendline), 0);
  if (recv(sockfd, recvline, MAXLINE,0) == 0){
   //error: server terminated prematurely
   perror("The server terminated prematurely"); 
   exit(4);
  }
  printf("%s", "String received from the server: ");
  fputs(recvline, stdout);
 }

#if (_WIN32) || (_WIN64)
 WSACleanup();
#endif

 exit(0);
}


Thursday, 14 May 2015

Neat and tidy FFT functions

Enjoy :


/* fft.c */

#include <math.h>


void fft_init (double *pFFTCoeffs,
    int FFTSize)

{
    register int i;

                /* Generate Sine and Cos tables */
    for (i = 0; i < ((3 * FFTSize) >> 2); i++)
    {
        *pFFTCoeffs++ = sin ((((double)i) * 6.283185307179586476925286766559) / ((double)FFTSize));
    }

}

void rfft (double RealData[],
    double ImagData[],
    const double *pFFTCoeffs,
    int FFTSize,
    int Log2Size)

{
    register const double   *pFFTSineCoeffs;
    register const double   *pFFTCosineCoeffs;
    register int            i, j, k, Stride, BflyCounter, g, h, HalfFFTSize;
    int                     Angle, AngleInc;    /* Angle step thro sin & cos tables */
    register double         RealTemp, ImagTemp, Cos, Sin;

    pFFTSineCoeffs = pFFTCoeffs;
    pFFTCosineCoeffs = pFFTCoeffs + (FFTSize >> 2);

    HalfFFTSize = FFTSize >> 1;

    Stride = j = HalfFFTSize;
    Angle = h = 0;

                        /* First stage */
    for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++)
    {
        RealTemp = RealData[h] - RealData[j];
        RealData[h] = RealData[h] + RealData[j];
        ImagData[h] = 0.0;      /* Clear imaginary part */
        RealData[j] = pFFTCosineCoeffs[Angle] * RealTemp;
        ImagData[j] = -pFFTSineCoeffs[Angle] * RealTemp;
        Angle++;
        h++;
        j++;
    }

    AngleInc = 2;

    for (i = 1; i < (Log2Size - 1); i++)    /* Middle stages */
    {
        k = Stride;
        Stride >>= 1;
        Angle = 0;
        for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++)
        {
            Cos = pFFTCosineCoeffs[Angle];
            Sin = pFFTSineCoeffs[Angle];
            Angle += AngleInc;

            h = BflyCounter;
            j = h + Stride;

            for (g = k; g <= FFTSize; g += k, h += k, j += k)
            {
                RealTemp = RealData[h] - RealData[j];
                ImagTemp = ImagData[h] - ImagData[j];
                RealData[h] = RealData[h] + RealData[j];
                ImagData[h] = ImagData[h] + ImagData[j];
                RealData[j] = Cos * RealTemp + Sin * ImagTemp;
                ImagData[j] = Cos * ImagTemp - Sin * RealTemp;
            }
        }
        AngleInc <<= 1;
    }

                        /* Final stage */
    for (h = 0, j = 1; h < FFTSize; h += 2, j += 2)
    {
        RealTemp = RealData[h] - RealData[j];
        ImagTemp = ImagData[h] - ImagData[j];
        RealData[h] = RealData[h] + RealData[j];
        ImagData[h] = ImagData[h] + ImagData[j];
        RealData[j] = RealTemp;     /* Cos = 1, sin = 0 */
        ImagData[j] = ImagTemp;
    }

                        /* Reorder scrambled data - slow mode uses less memory */
    for (j = 0, i = 0; j < FFTSize; j++)
    {
        if (j < i)
        {
            RealTemp = RealData[i]; RealData[i] = RealData[j]; RealData[j] = RealTemp;
            ImagTemp = ImagData[i]; ImagData[i] = ImagData[j]; ImagData[j] = ImagTemp;
        }

        k = HalfFFTSize;
        while ((k <= i) && (k >= 1))
        {
            i -= k;
            k >>= 1;
        }
        i += k;
    }

}   /* End of SDA_Rfft() */ 



void cfft (double RealData[],
    double ImagData[],
    const double *pFFTCoeffs,
    int FFTSize,
    int Log2Size)

{
    register int          i, j, k, Stride, BflyCounter, g, h;
    int                   Angle, AngleInc;    /* Angle step thro sin & cos tables */
    register double       ImagTemp, RealTemp, Cos, Sin;
    register const double *pFFTSineCoeffs = pFFTCoeffs;
    register const double *pFFTCosineCoeffs = pFFTCoeffs + (FFTSize >> 2);

    for (i = 0, Stride = FFTSize, AngleInc = 1; i < Log2Size; i++)
    {
        k = Stride;
        Stride >>= 1;
        Angle = 0;
        for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++)
        {
            Cos = pFFTCosineCoeffs[Angle];
            Sin = pFFTSineCoeffs[Angle];
            Angle += AngleInc;

            h = BflyCounter;
            j = h + Stride;

            for (g = k; g <= FFTSize; g += k, h += k, j += k)
            {
                RealTemp = RealData[h] - RealData[j];
                ImagTemp = ImagData[h] - ImagData[j];
                RealData[h] = RealData[h] + RealData[j];
                ImagData[h] = ImagData[h] + ImagData[j];
                RealData[j] = Cos * RealTemp + Sin * ImagTemp;
                ImagData[j] = Cos * ImagTemp - Sin * RealTemp;
            }
        }
        AngleInc <<= 1;
    }
    
                        /* Reorder scrambled data - slow mode uses less memory */
    for (j = 0, i = 0; j < FFTSize; j++)
    {
        if (j < i)
        {
            RealTemp = RealData[i]; RealData[i] = RealData[j]; RealData[j] = RealTemp;
            ImagTemp = ImagData[i]; ImagData[i] = ImagData[j]; ImagData[j] = ImagTemp;
        }

        k = FFTSize >> 1;
        while ((k <= i) && (k >= 1))
        {
            i -= k;
            k >>= 1;
        }
        i += k;
    }
}

These functions, and more, are available in the SigLib DSP Library.

Saturday, 28 February 2015

Debugging DSP Algorithms And Real-time Applications

I am often asked how I debug real-time and DSP functions. Here is a summary of what works for me.

While there are solutions for displaying signal data graphically, which we will discuss later, I often find that I need to log precise values so that I can check for single bit errors etc. To this end, printf is a pretty good solution but I find that a better one is to use fprintf to log the data to files, this way the tests can be automated and it is easy to use file comparison tools, like KDiff3, or python to check for differences between successive runs.

At a top level I’ll always use a
#define DEBUG_ON 1 // Set to 1 to enable debugging
and then use conditional compilation.

Next, I will use fprintf to send the results to a log file, called debug.log.
So the first thing we need to do is clear any existing open debug.log file, as follows :

int clearDebugfprintf (void)

{
    FILE *fp_LogFile;
    fp_LogFile = fopen( "debug.log", "w");
    if (fp_LogFile == NULL)
    {
        return (1);                  // Return error
    }
    fclose (fp_LogFile);

    return (0);                      // Return no error
}

For the actual logging I will typically use the following two functions or variations of them :

This function works in the same way as printf :

int debugfprintf (const char *ArgumentType, ...)

{
    va_list p_ArgumentList;

   FILE *fp_LogFile;
    fp_LogFile = fopen ("debug.log", "a");
    if (fp_LogFile == NULL)
    {
        return (1);
    }

    va_start (p_ArgumentList, ArgumentType);
    vfprintf (fp_LogFile, ArgumentType, p_ArgumentList);
    va_end (p_ArgumentList);
    fclose (fp_LogFile);

    return (0);
}


This function works in the same way as vprintf :

int debugvfprintf (const char *format,
    va_list ap)

{
    FILE *fp_LogFile;
    fp_LogFile = fopen( "debug.log", "a");
    if (fp_LogFile == NULL)
    {
        return (1);
    }

    vfprintf (fp_LogFile, format, ap);

    fclose (fp_LogFile);

    return (0);
}

One of the most common problems with these functions is how to include them in real-time applications without breaking the response time requirements. JTAG is horrendously slow for this requirement but on my favourite XMOS microcontrollers the debug data can be sent back to the development workstation using a parallel data path, called xSCOPE, in real-time - https://www.xmos.com/products/tools/xscope. xSCOPE uses the xCONNECT data path to route the data with very little overhead on the application being debugged.

To display data graphically I will typically use Gnuplot, as I documented in this blog post http://realgonegeek.blogspot.co.uk/2014/01/interfacing-cc-to-gnuplot.html, or the XMOS xSCOPE graphical display in xTIMEcomposer.

Saturday, 17 January 2015

DSP Tech Brief : How To Perform An FFT On A Real-Time Data Stream

In this example we are going to use the SigLib DSP Library functions to perform a Fast Fourier Transform on the Analog inputs on our laptop or desktop computer.

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

My laptop only supports a single microphone input to the ADC so to test the code I used an XMOS Multi-Function Audio (MFA) platform, shown in this image :



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 :-)

How To Connect Digital Signal Processing Functions To An ADC And A DAC

In my earlier post I showed how to implement an Infinite Impulse Response (IIR) filter :
Let's Do Some More Digital Signal Processing - Infinite Impulse Response (IIR) Filters.

In this blog I will look at how to interface DSP functions (e.g. the SigLib DSP library functions) to an Analog to Digital Converter (ADC) or Digital to Analog Converter (DAC). For this blog I will use the IIR filter functions mentioned in the above blog entry.

The system configuration is shown in this diagram :



In order to handle the ADC/DAC data input and output (I/O) it is necessary to write an Interrupt Service Routine (ISR) that responds to the interrupt from the ADC/DAC.

Any computer with a soundcard can be used to demonstrate this however modern Operating Systems trap interrupts and use drivers to interface to the hardware. In order to demonstrate hardware interrupts it is necessary to write a virtual software interrupt generator using a call back function.
I am going to use the free, cross-platform open-source audio I/O library PortAudio to implement the I/O then I will call the callback function that we will use as our ISR.

On a typical embedded system it will not be necessary to use PortAudio, just hook your ISR function to the interrupt from your ADC/DAC and it should work in exactly the same way.

Here is the header file (analog_io.h.) Note that I have declared the input and output registers as global variables. I know this doesn't conform to any standards for structured programming however these are designed to represent the input and output registers of the ADC and DAC devices.

int analog_open (int SampleRate, void (*analog_isr) (void));
int analog_close (void);

volatile short adc_in0, adc_in1, dac_out0, dac_out1;

Here is the wrapper code (analog_io.c.) for simulating our ISR using PortAudio. Note the pointer to the ISR is passed to the function analog_open () along with the desired sample rate.


#include <stdio.h>
#include "portaudio.h"
#include "analog_io.h"

#define PA_SAMPLE_TYPE      paInt32
#define FRAMES_PER_BUFFER   (64)

void (*p_analog_isr) (void);

typedef int SAMPLE;
PaStream *stream;

static int paCallback( const void *inputBuffer, void *outputBuffer,
                       unsigned long framesPerBuffer,
                       const PaStreamCallbackTimeInfo* timeInfo,
                       PaStreamCallbackFlags statusFlags,
                       void *userData );

/* PortAudio callback */
static int paCallback( const void *inputBuffer, void *outputBuffer,
                       unsigned long framesPerBuffer,
                       const PaStreamCallbackTimeInfo* timeInfo,
                       PaStreamCallbackFlags statusFlags,
                       void *userData )
{
    SAMPLE *out = (SAMPLE*)outputBuffer;
    const SAMPLE *in = (const SAMPLE*)inputBuffer;
    unsigned int i;
    (void) timeInfo; /* Prevent unused variable warnings. */
    (void) statusFlags;
    (void) userData;

    if( inputBuffer == NULL )         // If input buffer empty then output silence
    {
        for( i=0; i<framesPerBuffer; i++ )
        {
            *out++ = 0;               /* Left */
            *out++ = 0;               /* Right */
        }
    }
    else                              // If input buffer NOT empty then process data
    {
        for( i=0; i<framesPerBuffer; i++ )
        {
          adc_in0 = (short)(*in++ >> 16);            /* Left */
          adc_in1 = (short)(*in++ >> 16);            /* Right */
          p_analog_isr ();
          *out++ = (int)(dac_out0 << 16);
          *out++ = (int)(dac_out1 << 16);
        }
    }
    return paContinue;
}


int analog_open (int SampleRate, void (*analog_isr) (void))
{
  p_analog_isr = analog_isr;

  PaStreamParameters inputParameters, outputParameters;
  PaError err;

  err = Pa_Initialize();
  if( err != paNoError ) goto error;

  inputParameters.device = Pa_GetDefaultInputDevice(); /* default input device */
  if (inputParameters.device == paNoDevice) {
    fprintf(stderr,"Error: No default input device.\n");
    goto error;
  }
  inputParameters.channelCount = 2;       /* stereo input */
  inputParameters.sampleFormat = PA_SAMPLE_TYPE;
  inputParameters.suggestedLatency = Pa_GetDeviceInfo( inputParameters.device )->defaultLowInputLatency;
  inputParameters.hostApiSpecificStreamInfo = NULL;

  outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */
  if (outputParameters.device == paNoDevice) {
    fprintf(stderr,"Error: No default output device.\n");
    goto error;
  }
  outputParameters.channelCount = 2;       /* stereo output */
  outputParameters.sampleFormat = PA_SAMPLE_TYPE;
  outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
  outputParameters.hostApiSpecificStreamInfo = NULL;

  err = Pa_OpenStream(&stream,
                      &inputParameters,
                      &outputParameters,
                      SampleRate,
                      FRAMES_PER_BUFFER,
                      0, /* paClipOff, */  /* we won't output out of range samples so don't bother clipping them */
                      paCallback,
                      NULL );
  if( err != paNoError ) goto error;

  err = Pa_StartStream( stream );
  if( err != paNoError ) goto error;

  return 0;

error:
  Pa_Terminate();
  fprintf( stderr, "An error occured while using the portaudio stream\n" );
  fprintf( stderr, "Error number: %d\n", err );
  fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
  return -1;
}

int analog_close (void)

{
  int err;

  err = Pa_CloseStream( stream );
  if( err != paNoError ) goto error;

  Pa_Terminate();
  return 0;

error:
    Pa_Terminate();
    fprintf( stderr, "An error occured while using the portaudio stream\n" );
    fprintf( stderr, "Error number: %d\n", err );
    fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) );
  return -1;
}

Here is the header file (iir.h.) for the IIR filtering functions :

void iir_init (double *pState,
  const int NumberOfStages);

double iir_filter (const double Source,
  double *pState,
  const double *pCoeffs,
  const int NumberOfStages);

And here are the IIR filtering functions (taken from the blog mentioned earlier)  (iir.c.) :

#include "iir.h"

void iir_init (double *pState,
  const int NumberOfStages)

{
  int  i;

      /* Initialise the filter state array to 0 */
  for (i = 0; i < (NumberOfStages * 2); i++)
  {
    *pState++ = 0.0;
  }

}   /* End of iir_init () */


double iir_filter (const double Source,
  double *pState,
  const double *pCoeffs,
  const int NumberOfStages)

{
  register double       FeedbackSumOfProducts;
  register int i;
  register double       TempInputData;

  TempInputData = Source;

  for (i = 0; i < NumberOfStages; i++)
  {
    FeedbackSumOfProducts = TempInputData - (*(pCoeffs+3) * *pState) -
        (*(pCoeffs+4) * *(pState+1));       /* Feedback */
    TempInputData = (*pCoeffs * FeedbackSumOfProducts) +
        (*(pCoeffs + 1) * *pState) +
        (*(pCoeffs + 2) * *(pState+1));     /* Feedforward and save result for next time round */
    
    *(pState+1) = *pState;                  /* Move delayed samples */
    *pState = FeedbackSumOfProducts;

    pState += 2;        /* Increment array pointers */
    pCoeffs += 5;
  }

  return (TempInputData);                     /* Save output */

}   /* End of iir_filter () */

Now all we need is the top level program (pa_filter.c) which includes the filter coefficients (generated using Digital Filter Plus http://www.numerix-dsp.com/dfplus/), the main () function to initialize the filters and the analog interface and we need the ISR to call the DSP functions when the interrupt occurs.

#include <stdio.h>
#include <math.h>
#include "analog_io.h"
#include "iir.h"

#define SAMPLE_RATE     44100
#define FILTER_STAGES   1           /* Single 2nd ordrer biquad */

// Fc = 1 kHz LPF
double LPFCoefficientArray [] =
{
3.43448050764769180000e-003, 6.86896101529538370000e-003, 3.43448050764769180000e-003,
-1.89857508921973460000e+000, 9.12313011250325380000e-001
};

// Fc = 1 kHz HPF
double HPFCoefficientArray [] =
{
9.32663817644617140000e-001, -1.86532763528923430000e+000, 9.32663817644617140000e-001,
-1.85191320077338670000e+000, 8.78746746701696060000e-001
};

double Chan0FilterState [2 * FILTER_STAGES], Chan1FilterState [2 * FILTER_STAGES];


void analog_isr (void)
{
  double analog_sample;

                            // Process channel 0
  analog_sample = (double)adc_in0;
  analog_sample = iir_filter (analog_sample, Chan0FilterState, LPFCoefficientArray, FILTER_STAGES);
  dac_out0 = (short)(analog_sample);

                            // Process channel 1
  analog_sample = (double)adc_in1;
  analog_sample = iir_filter (analog_sample, Chan1FilterState, HPFCoefficientArray, FILTER_STAGES);
  dac_out1 = (short)(analog_sample);

}

int main(void)
{
  int Error;

  iir_init (Chan0FilterState, FILTER_STAGES);
  iir_init (Chan1FilterState, FILTER_STAGES);

  Error = analog_open (SAMPLE_RATE, analog_isr);  // Open the analog interface
  if (Error == -1)
    return 1;

  printf("Hit ENTER to stop program.\n");
  getchar();

  Error = analog_close ();                        // Close the analog interface
  if (Error == -1)
    return 1;

  return 0;
}

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_filter.c analog_io.c iir.c -DPA_USE_ASIO=1 portaudio_x86.lib

My laptop only supports a single microphone input to the ADC so to test the code I used an XMOS Multi-Function Audio (MFA) platform, shown in this image :



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 hear the difference between the two channels after applying the Low-pass and High-pass filters.