Sunday 21 January 2018

What's All This FFT Scaling Stuff Anyhow ?

In these earlier blog posts I briefly discussed FFT scaling however in this post I will look into it in more detail.



If I run the following Matlab or Octave code :

x = [0 0.7071 1 0.7071 0 -0.7071 -1 -0.7071];
y = fft(x)
z = ifft(y)

We see that y contains the expected DFT/FFT output (I.E. scaled by N) but z contains the iFFT output scaled to match the original time domain input :

y =
   0.0000 + 0.0000i   0.0000 - 4.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 4.0000i

z =
   0.00000   0.70710   1.00000   0.70710   0.00000  -0.70710  -1.00000  -0.70710

So Matlab and Octave both include additional operations on the output to multiply the result of the iFFT function by a factor of 1/N so that x == ifft(fft(x)) this makes a lot of sense and is documented in many DSP text books.

If we do the same thing in the SigLib DSP library :

Pseudo code :
    SLData_t x [] = {0., 0.7071, 1.,  0.7071,  0., -0.7071, -1., -0.7071};
    y = SDA_Rfft(x);
    z = SDA_Cifft(y);

We see the following output From SigLib :

y :
[0] =  0.000000 + j 0.000000
[1] =  0.000000 - j 3.999981
[2] =  0.000000 + j 0.000000
[3] =  0.000000 + j 0.000019
[4] =  0.000000 + j 0.000000
[5] =  0.000000 - j 0.000019
[6] =  0.000000 + j 0.000000
[7] = -0.000000 + j 3.999981

z :
[0] = 0.000000
[1] = 5.656800
[2] = 8.000000
[3] = 5.656800
[4] = 0.000000
[5] = -5.656800
[6] = -8.000000
[7] = -5.656800

This shows that the output of the iFFT has not been multiplied by 1/N. Now while his results in a different solution to Matlab, it also makes sense for real-time systems.
In a real-time system there are two considerations :

  1. In floating point systems, scaling can be left until the end of all of the processing and applied to give the correct result for the applications
  2. In fixed point systems, scaling is fundamental to achieving the best signal to noise performance for the algorithm, without overflowing the fixed point number system

If you do want to scale the SigLib output to match that of Matlab and Octave it is as simple as :

    SDA_Multiply (z,                        // Pointer to source array
                  SIGLIB_ONE / FFT_LENGTH,  // Inverse FFT length
                  z,                        // Pointer to destination array
                  FFT_LENGTH);              // Dataset length

And, of course, if this is a fixed point number system you can just right shift by log2(N) to achieve the same, which is very efficient.

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.

Evaluate The Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/

Saturday 20 January 2018

Compiling PortAudio with Cygwin/GCC

In this post : http://blog.numerix-dsp.com/2018/01/notes-on-using-cygwin-for-dsp.html. I described using Cygwin and GCC for developing Digital Signal Processing algorithms.

I like to use PortAudio for audio I/O and testing audio signal processing algorithms, particularly using the XMOS xCORE-200 Multichannel Audio Platform : https://www.xmos.com/support/boards?product=18334.

Here are my list of instructions for compiling and using PortAudio with Cygwin/GCC.

To build PortAudio using Cygwin/GCC, execute the following commands in the Cygwin command prompt (Note, the default configure/make creates a .dll library. This works perfectly well but I prefer static linking because then I do not have to worry about locating the .dll file) :

cd /cygdrive/c/portaudio_v19
chmod +x configure
./configure --enable-static --disable-shared --with-winapi=wmme,wasapi,wdmks
make

To get a list of the supported APIs, use the following command :

./configure --help

The results are :

--with-winapi           Select Windows API support
                          ([wmme|directx|asio|wasapi|wdmks][,...]) [wmme]

Now to build an example in a DOS or Cygwin command prompt (e.g. C:\portaudio_v19\examples\paex_sine.c) use the following command line :
gcc paex_sine.c -Wall -Wno-main -Wno-unused-value -std=c99 -I ../include -L ../lib/.libs -l portaudio -o paex_sine.exe -lrt -lwinmm -mthreads

pa_devs is a really useful utility for listing all the audio interfaces :
gcc pa_devs.c -Wall -Wno-main -Wno-unused-value -std=c99 -I ../include -L ../lib/.libs -l portaudio -o pa_devs.exe -lrt -lwinmm -mthreads

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.

Evaluate The Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/

Notes On Using Cygwin For DSP Development

Many companies I work with have corporate policies regarding the use of Windows OS so I usually develop application code in this OS. I have used a variety of compiler tools over the years but I have been increasingly using GCC under Cygwin for my DSP development work.

Over the past few years I have found the Cygwin/GCC/SigLb combination to be both powerful and reliable.

In order to assist customers setup a Cygwin/GCC/SigLb development environment I get the customer to install the base cygwin installation then run the following batch file :

rem Batch file for configuring Cygwin/GCC/SigLb development environment

rem Modify these environment variables

set CYGWIN_ROOT_TEMP=C:\cygwin64
set SIGLIB_ROOT_TEMP=/cygdrive/c/siglib

rem Do not modify below this line

setx CYGWIN_HOME %CYGWIN_ROOT_TEMP%
setx PATH "%PATH%;%CYGWIN_ROOT_TEMP%\bin"
setx PATH "%PATH%;C:\Program Files\gnuplot\bin;C:\siglib\Examples\FileIO"

setx C_INCLUDE_PATH %SIGLIB_ROOT_TEMP%/include:%SIGLIB_ROOT_TEMP%/nhl:%SIGLIB_ROOT_TEMP%/gnuplot_c/src
setx LIBRARY_PATH %SIGLIB_ROOT_TEMP%/lib/GCC_Win:%SIGLIB_ROOT_TEMP%/nhl:%SIGLIB_ROOT_TEMP%/gnuplot_c/src

rem Install required some packages (and some extras that are useful for development)

%CYGWIN_ROOT_TEMP%\setup-x86_64.exe --no-desktop --no-shortcuts --no-startmenu --quiet-mode --root %CYGWIN_ROOT_TEMP% --packages gcc-core gcc-g++ gccmakedep gdb git gitk gedit gedit-code-assistance gedit-devel gedit-plugins grep make cmake nano splint wget xorg-server xinit grep openssh diff patch

I use the Gnuplot/C (https://sourceforge.net/projects/gnuplotc/) interface to Gnuplot for plotting DSP graphical results. Gnuplot/C pipes the data from the application to Gnuplot and while Gnuplot can be installed in Cygwin I find that a native Windows installation is more efficient and there are no limitations when piping data from the Cygwin environment to the native Windows. More details on installing Gnuplot are available in this blog post : http://blog.numerix-dsp.com/2014/06/gnuplot-installation-under-windows.html.

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.

Evaluate The Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/

Wednesday 10 January 2018

Should I MAC Or Should I MSUB (With Apologies To The Clash)

When implementing FIR filters the coefficients are all in a linear array so the task of importing the coefficients into an application is relatively straight forward.

Importing IIR coefficients into an application, however, is quite different and there are several traps than can cause the application to fail.

The following diagram shows the flow diagram of a direct form I IIR biquad :


All digital filter design programs will return a single or multi-dimensional array along the lines of the following :

int pIIRCoeffs[NUMBER_OF_FILTER_STAGES * NUMBER_OF_COEFFS_PER_BIQUAD] = {
    Coefficients for Biquad #1,
    Coefficients for Biquad #2,
    Coefficients for Biquad #3
};

There are two traps to consider when using the coefficients from a digital filter design program.

  1. The sign of the feedback coefficients may or may not be negated, in which case you may need to negate the coefficients yourself.
  2. There is no international standard for which set of coefficients should be notated an or bn. You should check the documentation of your design tools because there are always : 3 feedforward coefficients and 2 feedback coefficients. You may need to re-order the feedforward and feedback coefficients in your array.
Point 1. raises further issues with the implementation :


  1. Negating the feedback coefficients allows the use of the Multiply-Accumulate (mac) instruction, which is available on all DSPs and many microprocessors and microcontrollers.
  2. Not negating the feedback coefficients will require a multiply followed by a subtract. Nearly all DSPs and some microprocessors and microcontrollers implement a Multiply-Subtract (msub) instruction. If your chosen device does not implement msub then the multiply and subtract operations will be separated leading to slower performance. In this case it would be better to negate the coefficients and use a code implementation that uses the mac instruction. Many digital filter design programs, including Numerix' Digital Filter Plus (http://www.numerix-dsp.com/dfplus/) provide an option to pre-negate the coefficients, which allows the use of the mac instruction. The SigLib DSP library (http://www.numerix-dsp.com/siglib.html) includes functions to implement either format.

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.

Evaluate The Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/

Tuesday 9 January 2018

Software For Editing Sampled Waveforms

I have talked before, in this blog, about using the excellent Audacity (http://audacity.sourceforge.net) for DSP applications.

When editing sampled data files and viewing them in the time or frequency domain then Audacity can be very powerful.

While Audacity does support some scripting functionality (https://manual.audacityteam.org/man/scripting.html), I often find that Sox (http://sox.sourceforge.net/) applications for pre-conditioning the audio files, extracting a section of audio or shortening a .wav file then I will typically use sox from a script.

On the other hand, if I want to do something more complex such as visually aligning the data in two files then Audacity is much more convenient for a small number of files.

For aligning a larger group of files I would typically do this using a correlator but that is a subject for a future blog post.

Evaluate the Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/.

New Version Of The SigLib DSP Library Released

Version 8.60 is the latest version of the SigLib DSP library and is available now from http://www.numerix-dsp.com/siglib.html.

Additional functionality includes new envelope detectors, new FIR and IIR filter architectures, microphone array beamformer and FFT functions.

To evaluate the Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/.