Monday, 23 December 2013

Arachnophobia !

This is just a little too freaky.

Dancing Robugtix T8 robot spider controlled by an XMOS xCORE device. :
http://www.youtube.com/watch?v=Yo2TUIEXJig

Enthusiastically reviewed by Adam Savage from Mythbusters :
http://www.youtube.com/watch?v=-vVblGlIMgw


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/

Sunday, 15 December 2013

How To Use MathGL Under Windows with the Microsoft Visual Studio 2010 64 bit compiler

Introduction

This solution is heavily based on the solution provided in the following tutorial : https://groups.google.com/forum/#!topic/mathgl/t-X9eg-5joU.
My thanks to RK for helping me find a solution.
I have found some changes in the packages since RK's post and I also needed 64 bit GLUT libraries so I have made minor some changes to his procedures.

Note : While this solution worked, I found that the MathGL functionality didn't quite stack up to that available in GNUPlot, see here for further details : http://realgonegeek.blogspot.co.uk/2014/01/interfacing-cc-to-gnuplot.html.

Requirements

  Microsoft Visual Studio 2010 : go.microsoft.com/fwlink/?LinkId=244366
  Microsoft SDK v7.1 64 bit compiler. Installation instructions are here : http://realgonegeek.blogspot.co.uk/2013/08/microsoft-visual-c-2010-sdk-v71-64-bit.html
  MathGL 2.2 source code : http://sourceforge.net/projects/mathgl/files/mathgl/mathgl%202.2/mathgl-2.2.tar.gz/download
  FreeGlut 2.8.1 (Martin Payne's Windows MSVC binaries : http://www.transmissionzero.co.uk/software/freeglut-devel/
  CMake 2.8.11 : http://www.cmake.org/files/v2.8/cmake-2.8.12.1-win32-x86.exe

Procedures
1. Follow steps 1 to 3 in RK's tutorial except using freeglut instead of the original GLUT.
2. Before running CMake-gui ensure that you have configured the Microsoft Compiler environment using :
  "C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\vcvarsall.bat" amd64
3. Follow steps 4 to 6 in RK's tutorial and choose "Visual Studio 10 Win64" in step 6.
4. Tick the Advanced box in CMake-gui to view all of the requried entries (e.g  GLUT_INCLUDE_DIR and GLUT_glut_LIBRARY)
5. Follow the rest of RK's tutorial
6. In your chosen source code folder create the following file (e.g. example.cpp) :

#pragma warning( disable: 4190 )    // Disable complex number warning

#include <mgl2/glut.h>

int sample(mglGraph *gr)
{
  mglData dat (100);
  for (long i = 0; i < 40; i++)
  {
    gr->NewFrame ();          // start frame
    gr->Box ();               // some plotting
    gr->Axis ();              // draw axis
    for (long j = 0; j < dat.nx; j++)
      dat.a[j] = sin (M_PI*j/dat.nx+M_PI*0.05*i);
    gr->Plot (dat, "b");      // "b" is colour ??
    gr->EndFrame ();          // end frame
    gr->WriteFrame ();        // save frame
  }
  return 0;
}

int main (int argc, char **argv)
{
  mglGLUT gr (sample, "MathGL examples");
  return 0;
}

7. Open a command window into your source folder and ensure the Microsoft Compiler environment is correctly configured, as per step 2. above.
8. Configure the compiler environment variables to locate the header and library files e.g. :
  set INCLUDE=%INCLUDE%;c:\mathgl-2.2\include
  set LIB=%LIB%;c:\mathgl_v2.2\src\Release;C:\mathgl_v2.2\widgets\Release
9. Compile example.cpp using the following command line :
  cl /EHsc example.cpp mgl-glut.lib
10. Copy the required dlls to the current folder :
  copy /Y C:\Users\johne2\Desktop\Graphics\mathgl\mathgl_v2.2\mathgl-2.2-cmake\src\Release\*.dll
  copy /Y C:\Users\johne2\Desktop\Graphics\mathgl\mathgl_v2.2\mathgl-2.2-cmake\widgets\Release\*.dll
  copy /Y C:\Users\johne2\Desktop\Graphics\mathgl\Glut\freeglut\bin\x64\*.dll
11. Execute example.exe and you should see a scrolling sine wave.
12. Play with MathGL :-)

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/

Friday, 13 September 2013

How To Access A Windows Library From Python, Using SWIG

In the previous blog (How To Build Windows Static And Dynamic Link Libraries) I described how to generate static and dynamic link libraries from C functions. This blog describes how to call these functions fom Python. In this example we will access the same VectorScalarMultiply function. We will use SWIG (http://www.swig.org) to provide the wrapper between Python and C. SWIG is a really useful tool for importing C libraries into almost all common scritping languages like Python, PERL, Java etc.

First off, we need to provide a SWIG interface file (Multiply.i) to tell SWIG about our library. SWIG is very powerful because it is able to obtain a lot of information from the library header file, for example, it can work out that Data_t is typedef'd to a double. However we do need to assist SWIG and tell it that we are going to be handling arrays of floating point data.
To do this we include the SWIG interface file "carrays.i" and tell SWIG we want to use the array_class interface. The complete interface file is shown here :

%module Multiply_wrap

%include "typemaps.i"
%include "carrays.i"

%array_class(double, doubleArray); /* Include support for "double *" functions */
%array_class(long, longArray); /* Include support for "long *" functions */

%{
#include "Multiply.h"
%}

%include "Multiply.h"

With this file we can now create the importable Python module that contains the functions we want to access. This is essentially a two stage process :
    Configure SWIG
    Compile the library into a Python importable module

Here is a complete batch file for creating a Python wrapper for accessing the library. Note that we are calling the new DLL (Multiply_wrap.dll) to differentiate it from Multiply.dll created in the earlier blog. Note that Python expects the file to be called Multiply_wrap.pyd so we will also rename it.
Please refer to the remarks for additional comments.

rem Delete the wrapper c file so that if it is not regenerated using swig then the compile will fail
rem this is helpful for development because if the you don't delete the old file then the compiler
rem will just process that, which might not be what you want !
del Multiply_wrap_wrap.c

c:\swig\swig -python -includeall Multiply_wrap.i

cl -Od -MD -DNDEBUG -DWIN32 -IC:\Anaconda\envs\py33\include -D_CONSOLE -DSTATIC_LIB=1 -DNO_STRICT Multiply_wrap_wrap.c /link user32.lib Multiply.lib C:\Anaconda\envs\py33\libs\python33.lib /DLL /out:Multiply_wrap.dll /NODEFAULTLIB:LIBCMT

rem The Python module should have extension .pyd
rem so we will delete the old one and rename the newly created dll
if exist _Multiply_wrap.pyd del /Q _Multiply_wrap.pyd

ren Multiply_wrap.dll _Multiply_wrap.pyd

Finally, we need to write a test script to access the C coded function.
This function includes some additional functionality that is not strictly necessary but this is to show some additional useful techniques such as :
    Printing the attributes of the library
    Converting Python lists to C arrays and back again
    Accessing the values of C #define constants from Python
Here is the complete Test.py script :

import multiply_wrap                        # Note the module name must be lower case !

print ('')
print ('Library attributes :')              # Print the attributes of the library
for a in dir(Multiply_wrap):
    print (a)
print ('')

a = [1.2, 3.4, 5.6, 7.8, 9.0]
b = [1000.1, 100.2, 1000.3, 1000.4, 1000.5] # Destination array to put results into

print ('')
print ('PI = ')                             # Print a constant defined in C
print (Multiply_wrap.PI)
print ('')

p_a = Multiply_wrap.doubleArray (5)
p_b = Multiply_wrap.doubleArray (5)

for index, item in enumerate(a):            # Copy data to array p_a
    p_a[index] = item

print ('a = ')                              # Print contents of a
for i in range(0,5):
    print(p_a[i])
print ('')

Multiply_wrap.VectorScalarMultiply (p_a, Multiply_wrap.TWO_PI, p_b, 5)

for i in range(0,5):                        # Copy data to b
    b[i] = p_b[i]

print ('Scaled a = ')
print (b)
print ('')

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/

Thursday, 12 September 2013

How To Build Windows Static And Dynamic Link Libraries With The Microsoft Compiler

Creating a C/C++ statically linked library is easy : compile the library code to object code (without creating an executable) then call the lib utility. In this case the same header file that declares the function can be used for both the library and the application.
Creating a dynamic linked library is a little bit more involved because the header file needs to be able to declare the functions for both "export" (from the DLL) and "import" (into the application). In this case the functions must be declared differently depending on whether they are being used in the library or the application.
Using some simple conditional compilation techniques allows the library functions to be written in a way that will support both static and dynamic linking.

In this example I am going to implement a common DSP function, vector scalar multiplication and call the function VectorScalarMultiply().

Using floating point data this function has the following declaration :

void VectorScalarMultiply (const double *,  /* Pointer to source array */
    const double,                           /* Multiplier */
    double *,                               /* Pointer to destination array */
    const long);                            /* Array length */

There are several things we can do to improve this function and the first is to make it portable to use different data types. We will use a typedef for this and the following code shows how to declare the typedef for double precision floating point numbers :

typedef double  Data_t;                 /* Data type */

Now, if we want to change the data type of the function all we need to do is change the typedef.

Secondly we can write this function declaration so that it is compilable into both static and dynamic libraries, as described above. For this we will use the following conditional compilation code :

#if defined (_MSC_VER)                  /* Defined by Microsoft compilers */
#include <windows.h>                    /* Required for Windows applications */

#if defined (STATIC_LIB)                /* Create statically linked library */
#define FUNC_DECL                       /* Not used with statically linked library */
#else                                   /* Create dynamically linked  library */
#if defined (DLL_SOURCE)                /* Defined on command line, if rebuilding DLL */
#define FUNC_DECL       __declspec(dllexport) WINAPI    /* DLL export function - used in DLL source */
#else
#define FUNC_DECL       __declspec(dllimport) WINAPI    /* DLL import function - used in Application */
#endif      // DLL_SOURCE
#endif      // STATIC_LIB
#endif      // _MSC_VER

In the above code, we will use a declaration on the compiler command line for the following values (-D is the Microsoft command line switch):
-DSTATIC_LIB=1
This is used to tell the compiler that the function is being compiled into a static library and the result is that the value FUNC_DECL becomes undeclared.
-DDLL_SOURCE=1
This is used to tell the compiler that the function is being compiled into a dynamic library and the result is that the value FUNC_DECL becomes "__declspec(dllexport) WINAPI".
The other option (if neither of these is declared) is : "__declspec(dllimport) WINAPI", which is used to include the header into the application.

This results in the following declaration :

void FUNC_DECL VectorScalarMultiply (const Data_t *,    /* Pointer to source array */
    const Data_t,                           /* Multiplier */
    Data_t *,                               /* Pointer to destination array */
    const long);                            /* Array length */

So the complete header file (which we will call Multiply.h) is as follows and we have defined a constant PI so that we can use it in the test application later on :

#define PI      3.141592

typedef double  Data_t;                 /* Data type */

#if defined (_MSC_VER)                  /* Defined by Microsoft compilers */
#include <windows.h>                    /* Required for Windows applications */

#if defined (STATIC_LIB)                /* Create statically linked library */
#define FUNC_DECL                       /* Not used with statically linked library */
#else                                   /* Create dynamically linked  library */
#if defined (DLL_SOURCE)                /* Defined on command line, if rebuilding DLL */
#define FUNC_DECL       __declspec(dllexport) WINAPI    /* DLL export function - used in DLL source */
#else
#define FUNC_DECL       __declspec(dllimport) WINAPI    /* DLL import function - used in Application */
#endif      // DLL_SOURCE
#endif      // STATIC_LIB
#endif      // _MSC_VER

#if defined (SWIG)                      /* Is this header included by SWIG ? */
#define FUNC_DECL
#endif

void FUNC_DECL VectorScalarMultiply (const Data_t *,    /* Pointer to source array */
    const Data_t,                           /* Multiplier */
    Data_t *,                               /* Pointer to destination array */
    const long);                            /* Array length */

Now we need to write the function, as follows, which we will put in a source file called Multiply.c :

#include "Multiply.h"

void FUNC_DECL VectorScalarMultiply (const Data_t * pSrc,
    const Data_t Multiplier,
    Data_t * pDst,
    const long SampleLength)

{
    register long   i;

    for (i = 0; i < SampleLength; i++)
    {
        *pDst++ = *pSrc++ * Multiplier;
    }

}       /* End of VectorScalarMultiply() */

This includes FUNC_DECL, for using the function in a static or dynamic link library. This also includes the data type declared as Data_t.

Finally we need to write a test application (Test.c), as follows :

#include "Multiply.h"
#include <stdio.h>

#define ARRAY_LENGTH    5L

void main (void)

{
    Data_t  a[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    Data_t  b[ARRAY_LENGTH];
    long     i;

    VectorScalarMultiply (a, PI, b, ARRAY_LENGTH);

    for (i = 0; i < ARRAY_LENGTH; i++)
    {
        printf ("b[%ld] = %lf\n", i, b[i]);
    }
}

Now all we need to do is compile and test our library and application. The appropriate command lines are as follows :

Static Linked Library :
To build the library :
cl /c -DSTATIC_LIB=1 Multiply.c
To add the module into the library (called Multiply.lib) :
lib Multiply.obj
To build the application :
cl -DSTATIC_LIB=1 Test.c Multiply.lib

Dynamic Linked Library :
To build the library :
cl /LD -DDLL_SOURCE=1 Multiply.c
This creates a dll called Multiply.dll and a library called Multiply.lib. The lib is a wrapper for the dll. It is entirely possible to access Multiply.dll directly from a C application but this is the easiest method.
To build the test application :
cl Test.c Multiply.lib

Now just execute Test.exe and the application should call the library function.

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/

Saturday, 31 August 2013

How To Plot Data From C Using Matplotlib

Following on from my earlier blog entitled "How To Pass A C Array To Python Solution" (http://realgonegeek.blogspot.co.uk/2013/08/how-to-pass-c-array-to-python-solution.html) here is some simple code to plot the same data using Matplotlib.

Note : While this solution worked, I found that it was very slow because of the time required to load Python so I ended up reverting back to GNUPlot, see here for further details : http://realgonegeek.blogspot.co.uk/2014/01/interfacing-cc-to-gnuplot.html.

All you need to do is take the earlier C example, change the Python method from PrintArray to PlotArray and recompile.

def PlotArray (a):
    import matplotlib.pyplot as pp
    from matplotlib.font_manager import FontProperties
    import numpy as np
    
    p1, = pp.plot (a, label = "line 1", color='red')

    ax = pp.subplot (111)
    box = ax.get_position ()
    ax.set_position ([box.x0, box.y0, box.width * 0.85, box.height]) # Shink x axis to fit legend
    fontP = FontProperties ()
    fontP.set_size ('small')

    pp.legend (loc=2, prop = fontP, fancybox = True, bbox_to_anchor=(1, 1.0))
    pp.title ("Title")
    pp.xlabel ("X-Axis")
    pp.ylabel ("Y-Axis")
    pp.ylim (ymin = 0, ymax = 10)
    pp.grid (True)

    pp.show (block=True)

    c = 0
    return c

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/

How To Pass A C Array To Python Solution

While I commonly use C or C++ for processing signals there are not too many options available for displaying the results graphically so I have been looking at using Matplotlib for the display. The first problem I encountered is how to pass an array of floating point data from C to Python.
There are several options and a lot of forum discussions on the subject but I couldn't find a simple solution. So here is one that I developed earlier.
This example works with an Anaconda/Python v3.3 installation with the 64 bit Microsoft compiler.

See here for an introduction to  Anaconda/Python v3.3 installation and calling Python from the 64 bit Microsoft compiler : http://realgonegeek.blogspot.co.uk/2013/08/how-to-embed-python-pylab-into-64bit.html.

Here is the C code :

#include <Python.h>
#include <stdio.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include "C:\Anaconda\envs\py33\Lib\site-packages\numpy\core\include\numpy\arrayobject.h"

double Array [] = {1.2, 3.4, 5.6, 7.8};

int main (int argc, char *argv[])
{
    PyObject *pName, *pModule, *pDict, *pFunc, *pArgs, *pValue;
    npy_intp dims[1] = { 4 };
    PyObject *py_array;

    Py_Initialize ();
    pName = PyUnicode_FromString ("PrintArray");
    // PyUnicode_FromString error checking here

    pModule = PyImport_Import (pName);                  // Load the module object
    // PyImport_Import error checking here
    Py_DECREF(pName);

    pFunc = PyObject_GetAttrString (pModule, "PrintArray");        // pFunc is a new reference
    // PyObject_GetAttrString error checking here

    import_array ();                                    // Required for the C-API : http://docs.scipy.org/doc/numpy/reference/c-api.array.html#importing-the-api

    py_array = PyArray_SimpleNewFromData (1, dims, NPY_DOUBLE, Array);
    // PyArray_SimpleNewFromData error checking here

    pDict = PyModule_GetDict (pModule);                 // pDict is a borrowed reference

    pArgs = PyTuple_New (1);
    PyTuple_SetItem (pArgs, 0, py_array);

    pFunc = PyDict_GetItemString (pDict, "PrintArray"); // pFunc is also a borrowed reference

    if (PyCallable_Check (pFunc)) 
    {
        PyObject_CallObject (pFunc, pArgs);
    } else 
    {
        printf ("Function not callable !\n");
    }

    Py_DECREF (py_array);                               // Clean up
    Py_DECREF (pModule);
    Py_DECREF (pDict);
    Py_DECREF (pFunc);

    Py_Finalize ();                                     // Finish the Python Interpreter

    return 0;
}

Here is the python script to print the contents of the array :

def PrintArray (a):
    print ("Contents of a :")
    print (a)
    c = 0
    return c

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/

How To Embed Python PyLab Into A 64bit Windows C Application

Here is a quick solution for embedding PyLab into a 64bit Windows C app. It was a bit fiddly to get setup so after a few trips down blind alleys, here is what I did.

The easiest way to get a Python installation that includes PyLab, Matplotlib, Numpy, Scipy etc is to download Anaconda from : http://www.continuum.io/downloads. In my case I wanted to use the Python v3.3 version under 64 bit Windows. The default Anaconda installation is for v2.x but it is very easy to update to Python v3.3 using conda. form C:\Anaconda type the following command :
conda create -n py33 python=3.3 anaconda

In order to use the v3.3 installation, prepend the following to the PATH environment variable :
C:\Anaconda\envs\py33;C:\Anaconda\envs\py33\Scripts;

Although Anaconda works very well from Python scripts I found that I need to set the following PYTHONPATH environment variable :
set PYTHONPATH=C:\Anaconda\envs\py33;C:\Anaconda\envs\py33\Lib;C:\Anaconda\envs\py33\DLLs

To configure the Visual C compiler for 64 bit mode :
call "C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\vcvarsall.bat" amd64
For further details on installing the 64 bit Microsoft compiler, please see here :
http://realgonegeek.blogspot.co.uk/2013/08/microsoft-visual-c-2010-sdk-v71-64-bit.html

To compile the app use the following command line :
cl example.c -IC:\Anaconda\envs\py33\include C:\Anaconda\envs\py33\libs\python33.lib /link /MACHINE:X64

If you want to try it, here is a quick test program :

#include "Python.h"
int main()
{
Py_Initialize ();
PyRun_SimpleString ("import sys");
PyRun_SimpleString ("import pylab");
PyRun_SimpleString ("print(sys.version)");
PyRun_SimpleString ("print()");
PyRun_SimpleString ("pylab.plot(range(10))");
PyRun_SimpleString ("pylab.show()");
Py_Exit(0);
return 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/

Tuesday, 27 August 2013

Microsoft Visual C++ 2010 / SDK v7.1 64 bit compiler installation failure - solution

I recently had to install the Microsoft Visual Studio 64 bit compiler on top of my existing Visual C++ 2010 installation, with the service pack SP1.
The 64 bit compiler is included in the Windows SDK 7.1.
If you try to install Windows SDK 7.1 on top of Visual C++ 2010 SP1 the installation process will fail.
Unfortunately the failure dialog that pops up after the installation failure gives no clue as to where the problem lies. You can open the failure log from the dialog but again the error messages do not provide a direct source of the problem.
Performing an internet search for solutions show up a lot of people having the same problem and a lot of other people not quite knowing the solution. One example of erroneous suggestion I came across was to install from the ISO rather than the web install.

Fortunately, those clever people at Mathworks have published a great set of instructions here to explain how to correctly install the 64 bit compiler :
http://www.mathworks.co.uk/support/solutions/en/data/1-ECUGQX/

A quick summary of the procedure is :
Uninstall the Visual C++ 2010 Redistributable packages
Install the SDK 7.1
Apply the compiler patch from Microsoft
Reinstall the Visual C++ 2010 Redistributable packages

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/

Monday, 1 July 2013

Digital Signal Processing - Bibliography

Here is a list of references that I have used over the years.

My goto reference for DSP is Oppenheim and Schafer's Discrete Time Signal Processing but for learning DSP I recommend Steven Smith's The Scientist and Engineer's Guide to Digital Signal Processing. Please see this blog entry :


As a point of historical note Park's Principles of Sigma-Delta Modulation for Analog-to-Digital Converters is a really useful reference for any narrow band signal processing applications and I have used it as a reference for wireless and ultrasound applications. If you can find a copy on the web it is well worth a read.

Analog Devices, DSP Applications Using The ADSP-2100 Family, Prentice Hall, USA.
Antoniou, A., (1993), Digital Filters : Analysis, Design and Applications, McGraw-Hill, USA.
Bellanger, M., (1984), Digital Processing Of Signals, Theory And Practice, John Wiley & Sons, UK
Burrus, C. S. and Parks, T. W., (1985), DFT/FFT and Convolution Algorithms, Wiley Interscience, NY, USA.
Cooley, J.W. and Tukey, J.W., (1965), An Algorithm for the Machine Calculation of Complex Fourier Series, Mathematics of Computation, The American Mathematical Society, USA.
Defatta, D. J., Lucas, J. G. and Hodgkiss, W. S., (1988), Digital Signal Processing - A system Design Approach, John Wiley & Sons Inc., USA.
Higgins, R. J., (1990), Digital Signal Processing in VLSI, Prentice-Hall, Englewood Cliffs, N.J., USA.
Harris, F.J., (1978), On The Use of Windows For Harmonic Analysis With the Discrete Fourier Transform, Proc. Inst. Electrical and Electronic Engineers, January.
Kay, Modern Spectral Estimation, Prentice-Hall Inc., USA
Knuth, D. E., (1981), Seminumerical Algorithms, vol 2 of The Art of Computer Programming, Addison-Wesley, USA.
Park S., Principles of Sigma-Delta Modulation for Analog-to-Digital Converters (APR8.PDF), Motorola Inc., USA
Marple, Digital Spectral Analysis with Applications, Prentice-Hall Inc., USA
Nussbaumer H., Fast Fourier Transform and Convolution Algorithms
Qian, S and Chen, D., (1996), Joint Time-Frequency Analysis, Prentice-Hall Inc., USA.
Oppenheim, A. V. and Schafer, R. W., (1989), Discrete Time Signal Processing, Prentice-Hall Inc., USA.
Orfanidis, Sophocles J., Introduction to Signal Processing, Prentice Hall, Inc. : http://eceweb1.rutgers.edu/~orfanidi/intro2sp/
Proakis, J. G., and Manolakis, D. G., (1992), Digital Signal Processing - Principles, Algorithms and Applications (2ed), Prentice-Hall, USA.
Proakis, J. G., Rader C. M., Ling, F. and Nikias, C. L., (1992), Advanced Digital Signal Processing, Macmillan Inc.
Smith, Steven W., The Scientist and Engineer's Guide to Digital Signal Processing : http://www.dspguide.com
Sorensen, H. V., Jones, D. L., Heideman, M. T. & Burrus, C. S., (1987), and Sig. Proc., Vol ASSP-35, No Real-valued Fast Fourier Transform Algorithms, IEEE Tran. Acoustics, Speech 6 pgs 849-863. (1981)
Stimson, G. W., (1998), Introduction To Airborne Radar (2ed), Scitech Publishing, USA.


Sunday, 30 June 2013

DSP Tech Brief : Some Notes On Implementing DSP Algorithms - C or Assembly, Fixed or Floating Point

There are many practical issues to consider when implementing DSP applications, including whether to choose fixed or floating-point devices, whether to program in a high level or assembler code. Here are some notes on the issues that may help guide a decision

One of the most popular techniques for developing DSP systems is to simulate the system in C on a general purpose micro-processor and then port the C code onto a DSP device, for many applications the C may give perfectly acceptable performance, in others however some of the routines must be coded in hand optimized assembler. Many of these functions can be easily optimized by compiling to assembler and then hand optimizing.

The choice between fixed and floating point devices is quite fundamental to an application. Fixed point devices generally have a higher clock speed due to their smaller size because they do not have the extra complexity of the floating point hardware. They are also produced in larger quantities which means that they benefit from the cost savings of mass production and they can therefore sell for a lower cost.

Floating-point devices give a higher dynamic range, which will allow more complex applications to be implemented for example a large 2D FFT can require a very large dynamic range, despite the limited resolution of the input, in this case it is almost certain that fixed point devices would not be suitable. One of the problems with fixed point devices is that the programmer must always be aware of the fact that the device can overflow its numerical bounds and if this is likely then the numbers should be scaled appropriately before continuing the process. The scaling can be quite time consuming and as a rule of thumb an FFT of more than about 1024 or 2048 points is quicker on a floating-point device, despite the faster clock speed of fixed point DSPs. Although floating-point devices give a larger dynamic range a floating-point device with a given mantissa length will give no greater precision than a fixed-point device with the same size word length.

When coding real-time routines in C there are some rules that can be followed to make the code as efficient as possible. First, all function parameters and local variables are placed on the stack and so they are accessed indirectly. This can often be over come by using compilers that use registers for variables but it is often the case that the code can be optimized by careful choice of register variables and by using global variables which are placed on the heap. When implementing interrupt service routines all registers that are used must be pushed onto the stack, to prevent side effects. In some cases using higher levels of optimization and hence extra registers for interrupt service routines may actually slow them down, due to the overhead of the extra stack manipulation that is required at both the start and end of the ISR.

Saturday, 29 June 2013

Another Good Way To Learn DSP – Online Courses

A while back I wrote an entry called : What Is The Best Way To Learn Digital Signal Processing. This listed a couple of signal processing related books that are available for download, for free. Another good way to learn DSP is to view some of the DSP Courses available on YouTube. In particular, these two are excellent. Don’t be put off by the fact that Alan’s is from 1975, fundamental theories don’t follow Moore’s Law ;-)


Of course, as good as these courses are, they don’t replace an instructor led course, like the ones from Oxford University ;-)

Friday, 8 March 2013

RaspBMC Hints And Tips

Having setup RaspBMC on my Pi I found that when playing music, that the display mode would go into fullscreen mode and I couldn't change the display mode using the standard TV remote control. With a little bit of playing I noticed that the Red, Green, Yellow and Blue buttons performed the following operations :
      Red : PVR

      Green : Videos
      Yellow : Music
      Blue : Pictures
I haven't got a PVR connected to the Pi so I figured I could reuse the Red button to swap in and out of fullscreen mode.

A quick web search showed that the key mappings are stored in the file :

/opt/xbmc-bcm/xbmc-bin/share/xbmc/system/keymaps/remote.xml
Loading this into nano and searching for "red" (without quotation marks) showed the following lines of code :

      <red>XBMC.ActivateWindow(MyPVR)</red>
      <green>XBMC.ActivateWindow(MyVideos)</green>
      <yellow>XBMC.ActivateWindow(MyMusic)</yellow>
      <blue>XBMC.ActivateWindow(MyPictures)</blue>
I then changed the <red> line to :
      <red>FullScreen</red>


I also wanted to change the fastforward and rewind buttons to skip forward and back in the tracklist so I made the following changes :
      <forward>FastForward</forward>
      <reverse>Rewind</reverse>
To :
      <forward>SkipNext</forward>
      <reverse>SkipPrevious</reverse>


I then saved the file and after a quick reboot everything worked as required.

I will update this post with new hints and tips as I find them.

Finally, 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.

Sunday, 24 February 2013

Raspberry Pi VPN Configuration


This blog entry explained how to setup a PPTPD VPN server on a Raspberry Pi and was based on a blog that is no longer available.

I recently had to install a VPN on another Pi so I thought it was about time to upgrade to the much more secure OpenVPN.

I followed the instructions here :
http://readwrite.com/2014/04/10/raspberry-pi-vpn-tutorial-server-secure-web-browsing
http://readwrite.com/2014/04/11/building-a-raspberry-pi-vpn-part-two-creating-an-encrypted-client-side
These instructions are excellent but I do advise that you proceed slowly because OpenVPN is tricky to setup.


Tuesday, 29 January 2013

How To Run A Command On A Remote Unix Machine From Windows

I recently had to create a script under Windows that could kick off a task under Unix. The user of the Windows machine had no knowledge of Unix so I wanted a script or application that I could put on the desktop of the windows machine, that the user could just double click to execute.

I decided to use Strawberry Perl (http://strawberryperl.com/) to execute a script under Windows. After installing Strawberry Perl I installed the perl module Net::SSH2 using the command :
$cpan Net::SSH2
Then I executed the following script :


#!/usr/local/bin/perl

#use strict;
#use warnings;
#use 5.010;  #if you are using version 5.10 or higher

use Net::SSH2;

my $host = "ipaddress";
my $user_id = "username";
my $password = "password";

my ($len, $buf); # Buffer for reading SSH2 output

my $ssh2 = Net::SSH2->new();
$ssh2->connect($host) or die "Unable to connect $host $@ \n";
$ssh2->auth_password($user_id,$password) or die "Unable to login $@ \n";
my $chan = $ssh2->channel();

#Initialize and execute Unix command remotely
$linuxCommand = "ls > output.tmp";
$chan->exec($linuxCommand) or die "Unable to execute the command $@ \n";

sleep (5);


In order to use this script in your application all you need to do is cut and paste the code above into an empty file with the extension .pl (e.g. remoteCommand.pl). Having done this, modify these lines :
my $host = "ipaddress";
my $user_id = "username";
my $password = "password";
(Do not remove the quotation marks)
Then modify this line to change the commande to execute on the remote machine :
$linuxCommand = "ls > output.tmp";
(Put your command between the quotation marks).

Now all you need to do is double click the Perl file and it will execute the command on the remote machine.

Finally, 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.

Friday, 18 January 2013

How to benchmark some C code (or assembly code, if you put the assembly in a C callable function)


I often have to benchmark DSP algorithms and while code profiling tools are excellent for large applications, DSP algorithms are often most easily benchmarked using simple timing functions.

Here are two examples that I interchange, depending on the underlying device and/or operating system. These examples are very portable but each solution has its own benefits and disadvantages and the two examples will almost always give different results when executed on the same platform/OS.

This first example uses the clock() function :


#include <stdio.h>
#include <time.h>

#define LOOP_COUNT 1000000000L

void main (void)
{
long i = 0L, j = 0L;

clock_t t_overhead, t_start, t_stop;

/* Calculate the timing overhead */
t_start = clock (); t_stop = clock (); t_overhead = t_stop - t_start;

t_start = clock ();
for (i = 0; i < LOOP_COUNT; i++) { j++; } // Put your code in here
t_stop = clock ();

printf ("Time : %lf seconds\n\n", (((double)(t_stop - t_start - t_overhead)) / CLOCKS_PER_SEC));

}


This example uses the ftime() function :


#include <stdio.h>
#include <sys/timeb.h>

#define LOOP_COUNT 1000000000L

void main (void)
{
long i = 0L, j = 0L;

struct timeb t_current;
long t_overhead, t_start, t_stop;

/* Calculate the timing overhead */
ftime(&t_current); t_start = (long)((t_current.time*1000L)+(long)t_current.millitm);
ftime(&t_current); t_stop = (long)((t_current.time*1000L)+(long)t_current.millitm);
t_overhead = t_stop - t_start;

ftime(&t_current); t_start = (long)((t_current.time*1000L)+(long)t_current.millitm);
for (i = 0; i < LOOP_COUNT; i++) { j++; } // Put your code in here
ftime(&t_current); t_stop = (long)((t_current.time*1000L)+(long)t_current.millitm);

printf ("Time : %lf seconds\n\n", (((double)(t_stop - t_start - t_overhead)) / 1000.));

}



Things to be aware of when using functions like this are :
While it is entirely possible to print times to the nearest ms I suspect that in many applications the values returned by the timing functions are not actually that accurate because underlying everything are a bunch of C functions, libraries and an OS.
When calling your processing function on a multi-tasking OS you have little or no control about what else the OS is doing in the background.
The clock function returns values that depend on the implementation
On some processors it returns 0 the first time it is called and then starts incrementing
On some processors it starts counting when the process starts
On some processors it counts in milliseconds
On some embedded processors (e.g. some DSPs I have programmed) clock () does not return the time in seconds but cycles of the processor.
You should always use CLOCKS_PER_SEC to turn this value into real time
The clock() function is more commonly implemented on embedded processors than the ftime() function.
If the function to be benchmarked is very short then it may be necessary to call it a number of times to average out the effects of the OS and the accuracy of the clock functions.



Friday, 11 January 2013

Raspberry Pi Hints And Tips


Actually, none of this is Pi specific but I have found that some friends have bought Pis but not had any previous experience with Linux so here are some notes I wrote for them.

To update the OS and applications and then clean up the install :
sudo apt-get update; sudo apt-get /y upgrade; sudo apt-get autoremove

To find out how much disk and memory space is used/available :
df -h; free -m

How to send an email from the command line :
echo -e "Subject: Test Mail\r\n\r\nThis is a test mail" |msmtp --debug --from=default -t destination@EmailAddress

How to add a new user :
sudo adduser username
sudo visudo
Add following line under the "root ALL=(ALL) ALL" Line:
username ALL=(ALL:ALL) NOPASSWD: ALL

You should delete the default user (user:raspberry, passwd:pi). How to delete a user :
Ensure you have backed up any important information then enter :
sudo deluser --remove-home username

How to change the hostname (e.g. from raspberrypi to donkey) :
sudo nano /etc/hostname
Change raspberrypi to donkey
sudo hostname donkey
nano /etc/hosts
change :
127.0.1.1       raspberrypi
to
127.0.1.1       donkey

Install Apache web server :
sudo apt-get install apache2 php5 libapache2-mod-php5

Tuesday, 8 January 2013

Introduction To Digital Signal Processing


Here is a quick introduction to DSP that I wrote many years ago, as part of a course : http://www.numerix-dsp.com/tutorials/DSP/DSPIntroduction.pdf.

Here is a quick introduction to frequency domain processing that I wrote around the same time : http://www.numerix-dsp.com/tutorials/DSP/FrequencyDomainProcessing.pdf.

Both of these are fairly high level so don't be disappointed if you are expecting the answer to life, the universe and everything.
That's 42, by the way ;-)

Sunday, 6 January 2013

Let's Do Some Digital Signal Processing (Part 2) - One Pole Filter


One-pole filters are very useful and efficient techniques for filtering signals. They use a feedback value to implement a low-pass filter. They can also be used to design efficient peak-hold filters in the time-domain and are also used in power spectrum analyzers and can help extract features of the signals that are being analyzed.

The structure of the one-pole filter is shown in the following diagram :
And the following equation :
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 :
Here is the code for the one-pole filter :

#include <stdio.h>

#define COUNT 8 /* Number of samples to process */
#define ALPHA 0.6 /* Feedback value */

float impulseSrc[] = {
 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 };
float stepSrc[] = {
 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
 };
float Dst[] = {
 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 };

void one_pole (const float *, float *, const float, float *, const int count);

void main (void)

{
 float State = 0.0; /* Reset the state value */
 one_pole (impulseSrc, Dst, 0.6, &State, COUNT);
 printf ("One Pole Impulse = %f, %f, %f, %f, %f, %f, %f, %f\n", Dst[0], Dst[1], Dst[2], Dst[3], Dst[4], Dst[5], Dst[6], Dst[7]);
 State = 0.0; /* Reset the state value */
 one_pole (stepSrc, Dst, 0.6, &State, COUNT);
 printf ("One Pole Step = %f, %f, %f, %f, %f, %f, %f, %f\n", Dst[0], Dst[1], Dst[2], Dst[3], Dst[4], Dst[5], Dst[6], Dst[7]);
}


void one_pole (const float *pSrc,
float *pDst,
const float Alpha,
float *pState,
const int VectorLength)
{
    float LocalState = *pState; /* Use a local variable to make the code more efficient */
    int   i;
    for (i = 0; i < VectorLength; i++)
    {
        LocalState = *pSrc++ + (LocalState * Alpha);
        *pDst++ = LocalState;
    }
    *pState = LocalState;
}

And here is the output :
One Pole Impulse = 1.000000, 0.600000, 0.360000, 0.216000, 0.129600, 0.077760, 0.046656, 0.027994
One Pole Step = 1.000000, 1.600000, 1.960000, 2.176000, 2.305600, 2.383360, 2.430016, 2.458010

This implementation shows one problem with feedback filters and that is the gain, as shown above with the step response. In an ideal scenario (and definitely when using fixed point devices) then the magnitude of the output should be the same as the magnitude of the input (i.e. 0 -> 1). In this example it can be seen that the output value (i.e. the gain) tends to 2.5 (i.e. 1/(1-Alpha)). To counter this we can rewrite the one-pole filter equation as follows :
y(n) = x(n) * (1- a) + a * y(n-1) 0 < a < 1.0

And the new flow diagram looks like :
Now when we rewrite the code, it looks like this :

void one_pole (float *pSrc,
float *pDst,
float Alpha,
         float *pState,
int VectorLength)
{
    float OneMinusAlpha;
    float LocalState = *pState;
    int   i;
    OneMinusAlpha = 1.0 - Alpha;

    for (i = 0; i < VectorLength; i++)
    {
        LocalState = ((LocalState * Alpha) + (*pSrc++ * OneMinusAlpha);
        *pDst++ = LocalState;
    }

    *pState = LocalState;
}

And now the gain is normalized.
As mentioned at the start of this blog, scaling is a big issue with feedback filters and the more feedback coefficients the bigger the problem. Solving the problem for one-pole filters is a trivial task, doing the same more biquad IIR filters is a much more complex task - maybe something for a future blog.