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.

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/

Copyright © 2015 Delta Numerix

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.

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/

Copyright © 2015 Delta Numerix

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.

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/

Copyright © 2015 Delta Numerix

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.

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/

Copyright © 2015 Delta Numerix

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

}

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/

Copyright © 2015 Delta Numerix

DSP Tech Brief : Neat and tidy FFT functions

Enjoy :


// fft.c

#include <math.h>


void fft_init (double *pFFTCoeffs,
    int FFTSize)

{
    int i;

                // Generate 3/4 Sine look-up-table
    for (i = 0; i < ((3 * FFTSize) >> 2); i++) {
        *pFFTCoeffs++ = sin ((((double)i) * 6.283185307179586476925286766559) / ((double)FFTSize));
    }

}   // End of fft_init()

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

{
    const double   *pFFTSineCoeffs;
    const double   *pFFTCosineCoeffs;
    int            i, j, k, Stride, BflyCounter, g, h, HalfFFTSize;
    int            Angle, AngleInc;    // Angle step thro sin & cos tables
    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
    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 rfft() 



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

{
    int          i, j, k, Stride, BflyCounter, g, h;
    int          Angle, AngleInc;    // Angle step thro sin & cos tables
    double       ImagTemp, RealTemp, Cos, Sin;
    const double *pFFTSineCoeffs = pFFTCoeffs;
    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
    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;
    }
}   // End of cfft()

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

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/

Copyright © 2015 Delta Numerix

DSP Tech Brief : 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.

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/

Copyright © 2015 Delta Numerix