Saturday, 28 July 2018

How To Interpolate The Output Of An FFT To Get Higher Resolution

One common requirement for a frequency analyzer is to maximize the amplitude or frequency estimation estimate for a given sinusoid.

One easy way to do this is to increase the size of the FFT however this comes at the cost of reduced time domain resolution and increased requirement for MIPS and memory.

A low impact solution to increase the resolution of the output is to interpolate the vertex to give a better estimate of the amplitude and frequency of an individual sinusoid.

The overall procedure is shown in the following diagram :

FFT Interpolation Functions

Following the FFT function is a peak detector followed by a peak interpolator.
The interpolation function uses 3 point vertex quadratic interpolation, as shown in the following diagram.

Peak Interpoation

Given that xm is the location of the peak and ym is the magnitude of the peak we can interpolate the peak value efficiently by setting xm-1 == 0, xm == 1, xm+1 == 2, which allows the interpolation to be implemented using the following equations :

a = (ym-1 - ym) * (ym-1 - ym+1)
b = -2 * (ym-1 -ym) + 0.5 * (ym-1 - ym+1)

Interpolated Frequency = -0.5 * b / a
Interpolated Magnitude = ym-1 - 0.25 * b2 / a

To analyze the results, we ran a series of sinusoids through the algorithm and plotted the percentage error. In the following graph we have a normalized sample rate of 1.0 Hz and we are using a 512 point FFT. Plotting the error percentage results between FFT bin 7 (0.013671875 Hz) and FFT bin 9 (0.017578125 Hz).

Interpolation Function Results

We can observe that the maximum gain error percentage = 4.519154% and the maximum frequency error percentage = 0.596755%.

Note, this solution is very accurate when the signal constitutes a single sinusoid or a sinusoid sampled with low noise. However the accuracy will vary with increasing noise an/or adjacent sinusoids - chosing a good windowing function will definitely help with the latter.

Functions to implement the FFT interpoator are included in the SigLib DSP library. To evaluate the Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/.

Friday, 22 June 2018

Really Useful Digital Filter Design Spreadsheet

I have uploaded a Really Useful Digital Filter Design Spreadsheet to my Applications Note section : http://www.numerix-dsp.com/appsnotes/ because it no longer appears on the Music-DSP archive, where I originally found it.

It implements the equations from the excellent RBJ Audio-EQ-Cookbook and has some minor fixes.

I have used this spreadsheet for a lot of audio filter design work and also teaching students. I have tried to track down the original author, to give them the credit they deserve. If anyone knows who wrote this spreadsheet then please do ask them to get in touch.

Tuesday, 20 March 2018

New Version Of The SigLib DSP Library Released

Version 8.61 is the latest version of the SigLib Digital Signal Processing (DSP) library and is available now from http://www.numerix-dsp.com/siglib.html.

V8.61 now includes support for Android OS.

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

Monday, 19 March 2018

How To Build A Shared Library Under Android Studio And Import Into A Project

I recently wanted to build the SigLib Digital Signal Processing library into an Android library and import into a java application.

My preferred API for cross-language integration is SWIG.

While there are many tutorials available to explain how to create an integrated Android Studio project that includes java and native C++, I struggled to find one that treated the two components separately.

Here is my solution.

Build The Library From The Sources

I used Window 10 as the main development environment and Android Studio for the main application building however I found that it was easier to build the library on the command line under Cygwin.

Install Android Studio from here : https://developer.android.com/studio/install.html

Install NDK, CMake and LLDB from within Android Studio using :
    File | Settings | SDK Tools
    Then tick the boxed and click OK

When you try to build a C program that includes header files such as stdio.h you will receive an error message along the lines of :

  /home/John/Android/Sdk/ndk-bundle/sysroot/usr/include/linux/types.h:21:10: fatal error: 'asm/types.h' file not found
  #include <asm/types.h>
           ^~~~~~~~~~~~~
  1 error generated.

The solution is here : https://github.com/android-ndk/ndk/issues/510
Modify setup-toolchain.mk as follows :

--- a/build/core/setup-toolchain.mk
+++ b/build/core/setup-toolchain.mk
@@ -152 +152 @@
-    -isystem $(SYSROOT_INC)/usr/include/$(header_triple_$(TARGET_ARCH))
+    -isystem $$(call host-path,$(SYSROOT_INC)/usr/include/$(header_triple_$(TARGET_ARCH)))

Open up a Cygwin command prompt and cd to the folder with the library source. Set the path environment variables to the tools :

export PATH=$PATH:/cygdrive/c/Users/John/AppData/Local/Android/Sdk/ndk-bundle/toolchains/arm-linux-androideabi-4.9/prebuilt/windows-x86_64
export PATH=$PATH:/cygdrive/c/Users/John/AppData/Local/Android/Sdk/ndk-bundle/build

Build an application or library :

ndk-build NDK_PROJECT_PATH=$(pwd) APP_BUILD_SCRIPT=$(pwd)/Android.mk

That should generate the .so and all the necessary .java wrapper files.

Adding The Library To An Android Studio Project

For my project, I started with the simple GraphView (http://www.android-graphview.org/) application : described in this video https://www.youtube.com/watch?v=zbTvJZX0UDk. The source code is available here : https://github.com/mitchtabian/CreateSimpleGraphView. Thanks Mitch for the great tutorial.


Added to app/build.gradle

    sourceSets.main {
        jni.srcDirs = [] // This prevents the auto generation of Android.mk
        jniLibs.srcDir 'siglib_wrap' // Import precompiled libraries in your project.
    }


copy full SWIG installation, from above, to : project/app/src/main/java/com/numerix_dsp/siglib_wrap

copy .so to : project/app/src/main/jniLibs


References

[1] The Main Reference : http://www.sureshjoshi.com/mobile/android-ndk-in-android-studio-with-swig/

[2] SWIG Android : http://orx-project.org/wiki/en/tutorials/community/enobayram/swig_android

[3] SWIG Wrapping : https://stackoverflow.com/questions/30110579/how-do-i-wrap-an-existing-c-library-for-use-with-android-studio-with-swig-give

[4] JNI : https://stackoverflow.com/questions/21096819/jni-and-gradle-in-android-studio/26693354#26693354

[5] SWIG : http://www.swig.org/Doc3.0/
    http://www.swig.org/Doc3.0/Android.html
    http://www.swig.org/Doc3.0/Java.html

[6] Android Libraries : https://developer.android.com/studio/projects/android-library.html

[7] Android C/C++ : https://developer.android.com/studio/projects/add-native-code.html

[8] Compiling Libraries For Android : https://tariqzubairy.wordpress.com/2012/03/09/arm-binaries-static-library-for-android/

Application Source Code

package com.numerix_dsp.siglib_graph;

import android.os.Bundle;

import android.support.v7.app.AppCompatActivity;

import com.jjoe64.graphview.GraphView;
import com.jjoe64.graphview.series.DataPoint;
import com.jjoe64.graphview.series.LineGraphSeries;

import siglib_wrap.siglib_wrap;
import siglib_wrap.siglib_wrapJNI;
import siglib_wrap.SLSignal_t;
import siglib_wrap.SLSignalFillMode_t;
import siglib_wrap.SWIGTYPE_p_double;

public class SigLib_Graph extends AppCompatActivity {

    private siglib_wrap mSigLib;

    static {
        try {
            System.loadLibrary("siglib_wrap");
        } catch (UnsatisfiedLinkError e) {
            System.err.println("siglib_wrap library failed to load.\n" + e);
            System.exit(1);
        }
    }

    int SAMPLE_LENGTH = 512;

    LineGraphSeries<DataPoint> series;

    @Override
    protected void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        setContentView(R.layout.activity_sig_lib__graph);

        mSigLib = new siglib_wrap();

        SWIGTYPE_p_double SinePhase = mSigLib.new_doubleArray(1);
        mSigLib.doubleArray_setitem(SinePhase, 0, 0.0);

        SWIGTYPE_p_double timeArray = mSigLib.new_doubleArray(SAMPLE_LENGTH);
        SWIGTYPE_p_double nullArray = mSigLib.new_doubleArray(1);

        mSigLib.SDA_SignalGenerate (timeArray,          // Output array pointer
                SLSignal_t.SIGLIB_SINE_WAVE,            // Signal type - Sine wave
                siglib_wrapJNI.SIGLIB_ONE_get(),        // Signal peak level
                SLSignalFillMode_t.SIGLIB_FILL,         // Fill (overwrite) or add to existing array contents
                0.015,                              // Signal frequency
                siglib_wrapJNI.SIGLIB_ZERO_get(),       // D.C. Offset
                siglib_wrapJNI.SIGLIB_ZERO_get(),       // Unused
                siglib_wrapJNI.SIGLIB_ZERO_get(),       // Signal end value - Unused
                SinePhase,                              // Signal phase - maintained across array boundaries
                nullArray,                              // Unused
                SAMPLE_LENGTH);                         // Output dataset length

        GraphView graph = (GraphView) findViewById(R.id.graph1);
        series = new LineGraphSeries<DataPoint>();
        for(int i =0; i<SAMPLE_LENGTH; i++) {
            series.appendData(new DataPoint((double)i, siglib_wrap.doubleArray_getitem(timeArray, i)), true, SAMPLE_LENGTH);
        }
        graph.addSeries(series);
    }
}

The full Android Studio project can be downloaded as part of the SigLib evaluation package, from here : http://www.numerix-dsp.com/eval/

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

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/