Blog

Archive for the ‘NMath’ Category

Filtering with Wavelet Transforms

Friday, December 18th, 2015

Discrete time wavelet transforms have found engineering applications in computer vision, pattern recognition, signal filtering and perhaps most widely in signal and image compression. In 2000 the ISO JPEG committee proposed a new JPEG2000 image compression standard that is based on the wavelet transform using two Daubechies wavelets. This standard made the relatively new image decomposition algorithm ubiquitous on desktops around the world.

In signal processing, wavelets have been widely investigated for use in filtering bio-electric signals, among many other applications. Bio-electric signals are good candidates for multi-resolution wavelet filtering over standard Fourier analysis due to their non-stationary character. In this article we’ll discuss the filtering of electrocardiograms or ECGs and demonstrate with code examples how to filter an ECG waveform using NMath‘s new wavelet classes; keeping in mind that the techniques and code shown here apply to a wide class of time series measurements. If wavelets and their applications to filtering are unfamiliar to the reader, read a gentle and brief introduction to the subject in Wavelets for Kids: A Tutorial Introduction [5].

Filtering a Time Series with Wavelets

PhysioNet provides free access to a large collections of recorded physiologic signals, including many ECG’s. The ECG signal we will filter here, named aami-ec13 on PhysioNet, is shown below.

ECG Signal

ECG Signal

Our goal will be to remove the high frequency noise while preserving the character of the wave form, including the high frequency transitions at the signal peaks. Fourier based filter methods are ill suited for filtering this type of signal due to both it’s non-stationarity, as mentioned, but also the need to preserve the peak locations (phase) and shape.

A Wavelet Filter

As with Fourier analysis there are three basic steps to filtering signals using wavelets.

  1. Decompose the signal using the DWT.
  2. Filter the signal in the wavelet space using thresholding.
  3. Invert the filtered signal to reconstruct the original, now filtered signal, using the inverse DWT.

Briefly, the filtering of signals using wavelets is based on the idea that as the DWT decomposes the signal into details and approximation parts, at some scale the details contain mostly the insignificant noise and can be removed or zeroed out using thresholding without affecting the signal. This idea is discussed in more detail in the introductory paper [5]. To implement this DWT filtering scheme there are two basic filter design parameters: the wavelet type and a threshold. Typically the shape and form of the signal to be filtered is qualitatively matched to the general shape of the wavelet. In this example we will use the Daubechies forth order wavelet.

ECG Signal closeup

ECG Waveform

Daubechies 4 wavelet

Daubechies 4 wavelet

The general shape of this wavelet roughly matches, at various scales, the morphology of the ECG signal. Currently NMath supports the following wavelet families: Harr, Daubechies, Symlet, Best Localized, and Coiflet, 27 in all. Additionally, any custom wavelet of your invention can be created by passing in the wavelet’s low & high pass decimation filter values. The wavelet class then imposes the wavelet’s symmetry properties to compute the reconstruction filters.

   // Build a Coiflet wavelet.
   var wavelet = new FloatWavelet( Wavelet.Wavelets.C4 );
 
   // Build a custom reverse bi-orthogonal wavelet.
   var wavelet = new DoubleWavelet( new double[] {0.0, 0.0, 0.7071068, 0.7071068, 0.0, 0.0}, new double[] {0.0883883, 0.0883883, -0.7071068, 0.7071068, -0.0883883, -0.0883883} );

The FloatDWT class provides four different thresholding strategies: Universal, UniversalMAD, Sure, and Hybrid (a.k.a SureShrink). We’ll use the Universal threshold strategy here. This is a good starting point but this strategy can over smooth the signal. Typically some empirical experimentation is done here to find the best threshold for the data (see [1], also see [4] for a good overview of common thresholding strategies.)

Wavelet Filtering Code

The three steps outlined above are easily coded using two classes in the NMath library: the FloatDWT class and the FloatWavelet class. As always in NMath, the library offers both a float precision and a double precision version of each of these classes. Let’s look at a code snippet that implements a DWT based filter with NMath.

   // Choose wavelet, the Daubechies 4 wavelet
   var wavelet = new FloatWavelet( Wavelet.Wavelets.D4 );
 
   // Build DWT object using our wavelet & data
   var dwt = new FloatDWT( data, wavelet );
 
   // Decompose signal with DWT to level 5
   dwt.Decompose( 5 );
 
   // Find Universal threshold & threshold all detail levels
   double lambdaU = dwt.ComputeThreshold( FloatDWT.ThresholdMethod.Universal, 1 );
   dwt.ThresholdAllLevels( FloatDWT.ThresholdPolicy.Soft, new double[] { lambdaU, 
       lambdaU, lambdaU, lambdaU, lambdaU } );
 
   // Rebuild the filtered signal.
   float[] reconstructedData = dwt.Reconstruct();

The first two lines of code build the wavelet object and the DWT object using both the input data signal and the abbreviated Daubechies wavelet name Wavelet.Wavelets.D4. The third line of code executes the wavelet decomposition at five consecutive scales. Both the signal’s details and approximations are stored in the DWT object at each step in the decomposition. Next, the Universal threshold is computed and the wavelet details are thresholded using the same threshold with a Soft policy (see [1], pg. 63). Lastly, the now filtered signal is reconstructed.

Below, the chart on the left shows the unfiltered ECG signal and the chart on the right shows the wavelet filtered ECG signal. It’s clear that this filter very effectively removed the noise while preserving the signal.

Raw ECG Signal

Raw ECG Signal

Filtered ECG Signal

Filtered ECG Signal

These two charts below show a detail from the chart above from indicies 500 to 1000. Note how well the signal shape, phase, and amplitude has been preserved in this non-stationary wavelet-filtered signal.

Detail Raw ECG Signal

Detail Raw ECG Signal

Detail Filtered ECG Signal

Detail Filtered ECG Signal

It is this ability to preserve phase, form, and amplitude in DWT based filters all while having a O(n log n) runtime that Fourier-based filters enjoy that has made wavelets such an important part of signal processing today. The complete code for this example along with a link to the ECG data is provided below.

Paul

References

[1] Guomin Luo and Daming Zhang (2012). Wavelet Denoising, Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology, Dr. Dumitru Baleanu (Ed.), ISBN: 978-953-51-0494-0, InTech, pp. 59-80. Available from: http://www.intechopen.com/books/advances-in-wavelet-theory-and-their-applicationsin-engineering-physics-and-technology/wavelet-denoising
[2] Burhan Ergen (2012). Signal and Image Denoising Using Wavelet Transform, Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology, Dr. Dumitru Baleanu (Ed.), ISBN: 978-953-51-0494-0, InTech, DOI: 10.5772/36434. Available from: http://www.intechopen.com/books/advances-in-wavelet-theory-and-their-applications-in-engineering-physics-and-technology/wavelet-signal-and-image-denoising
[3] Rami Cohen: Signal Denoising Using Wavelets, Project Report, 2012. Available from: http://tx.technion.ac.il/~rc/SignalDenoisingUsingWavelets_RamiCohen.pdf
[4] M. C. E. Rosas-Orea, M. Hernandez-Diaz, V. Alarcon-Aquino, and L. G. Guerrero-Ojeda, A Comparative Simulation Study of Wavelet Based Denoising Algorithms, Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005), 2005 © IEEE
[5] Brani Vidakovic and Peter Mueller, Wavelets for Kids: A Tutorial Introduction, Duke University, 1991. Available from: http://gtwavelet.bme.gatech.edu/wp/kidsA.pdf

Test Data

To copy the data file provided by PhysioNet for this example click: ECG_AAMIEC13.data
This ECG data was taken from the ANSI EC13 test data set waveforms.

Complete Test Code

 
    public void BlogECGExample()
    {
      // Define your own dataDir
      var dataDir = "................";
 
      // Load ECG wave from physionet.org data file.
      string filename = Path.Combine( dataDir, "ECG_AAMIEC13.data.txt" );
      string line;
      int cnt = 0;
      FloatVector ecgMeasurement = new FloatVector( 3000 );
      var fileStrm = new System.IO.StreamReader( filename );
      fileStrm.ReadLine(); fileStrm.ReadLine();
      while ( ( line = fileStrm.ReadLine() ) != null && cnt < 3000 )
      {
        ecgMeasurement[cnt] = Single.Parse( line.Split( ',' )[1] );
        cnt++;
      }
 
      // Choose wavelet
      var wavelet = new FloatWavelet( Wavelet.Wavelets.D4 );
 
      // Build DWT object
      var dwt = new FloatDWT( ecgMeasurement.DataBlock.Data, wavelet );
 
      // Decompose signal with DWT to level 5
      dwt.Decompose( 5 );
 
      // Find Universal threshold & threshold all detail levels with lambdaU
      double lambdaU = dwt.ComputeThreshold( FloatDWT.ThresholdMethod.Universal, 1 );
      dwt.ThresholdAllLevels( FloatDWT.ThresholdPolicy.Soft, new double[] { lambdaU, lambdaU, lambdaU, lambdaU, lambdaU } );
 
      // Rebuild the signal to level 1 - the original (filtered) signal.
      float[] reconstructedData = dwt.Reconstruct();
 
      // Display DWT results.
      BlogECGExampleBuildCharts( dwt, ecgMeasurement, reconstructedData );
 
    }
 
    public void BlogECGExampleBuildCharts( FloatDWT dwt, FloatVector ECGMeasurement, float[] ReconstructedData )
    {
 
      // Plot out approximations at various levels of decomposition.
      var approxAllLevels = new FloatVector();
      for ( int n = 5; n > 0; n-- )
      {
        var approx = new FloatVector( dwt.WaveletCoefficients( DiscreteWaveletTransform.WaveletCoefficientType.Approximation, n ) );
        approxAllLevels.Append( new FloatVector( approx ) );
      }
 
      var detailsAllLevels = new FloatVector();
      for ( int n = 5; n > 0; n-- )
      {
        var approx = new FloatVector( dwt.WaveletCoefficients( DiscreteWaveletTransform.WaveletCoefficientType.Details, n ) );
        detailsAllLevels.Append( new FloatVector( approx ) );
      }
 
      // Create and display charts.
      Chart chart0 = NMathChart.ToChart( detailsAllLevels );
      chart0.Titles.Add( "Concatenated DWT Details to Level 5" );
      chart0.ChartAreas[0].AxisY.Title = "DWT Details";
      chart0.Height = 270;
      NMathChart.Show( chart0 );
 
      Chart chart1 = NMathChart.ToChart( approxAllLevels );
      chart1.Titles.Add("Concatenated DWT Approximations to Level 5");
      chart1.ChartAreas[0].AxisY.Title = "DWT Approximations";
      chart1.Height = 270;
      NMathChart.Show( chart1 );
 
      Chart chart2 = NMathChart.ToChart( (new FloatVector( ReconstructedData ))[new Slice(500,500)] );
      chart2.Titles[0].Text = "Thresholded & Reconstructed ECG Signal";
      chart2.ChartAreas[0].AxisY.Title = "mV";
      chart2.Height= 270;
      NMathChart.Show( chart2 );
 
      Chart chart3 = NMathChart.ToChart( (new FloatVector( ECGMeasurement ))[new Slice(500,500)] );
      chart3.Titles[0].Text = "Raw ECG Signal";
      chart3.ChartAreas[0].AxisY.Title = "mV";
      chart3.Height = 270;
      NMathChart.Show( chart3 );
 
    }

Precision and Reproducibility in Computing

Monday, November 16th, 2015

Run-to-run reproducibility in computing is often assumed as an obvious truth. However software running on modern computer architectures, among many other processes, particularly when coupled with advanced performance-optimized libraries, is often only guaranteed to produce reproducible results only up to a certain precision; beyond that results can and do vary run-to-run. Reproducibility is interrelated with the precision of floating-point point types and the resultant rounding, operation re-ordering, memory structure and use, and finally how real numbers are represented internally in a computer’s registers.

This issue of reproducibility arises with NMath users when writing and running unit tests; which is why it’s important when writing tests to compare floating point numbers only up to their designed precision, at an absolute maximum. With the IEEE 754 floating point representation which virtually all modern computers adhere to, the single precision float type uses 32 bits or 4 bytes and offers 24 bits of precision or about 7 decimal digits. While the double precision double type requires 64 bits or 8 bytes and offers 53 bits of precision or about 15 decimal digits. Few algorithms can achieve significant results to the 15th decimal place due to rounding, loss of precision due to subtraction and other sources of numerical precision degradation. NMath’s numerical results are tested, at a maximum, to the 14th decimal place.

A Precision Example

As an example, what does the following code output?

      double x = .050000000000000003;
      double y = .050000000000000000;
      if ( x == y )
        Console.WriteLine( "x is y" );
      else
        Console.WriteLine( "x is not y" );

I get “x is y”, which is clearly not the case, but the number x specified is beyond the precision of a double type.

Due to these limits on decimal number representation and the resulting rounding, the numerical results of some operations can be affected by the associative reordering of operations. For example, in some cases a*x + a*z may not equal a*(x + z) with floating point types. Although this can be difficult to test using modern optimizing compilers because the code you write and the code that runs can be organized in a very different way, but is mathematically equivalent if not numerically.

So reproducibility is impacted by precision via dynamic operation reorderings in the ALU and additionally by run-time processor dispatching, data-array alignment, and variation in thread number among other factors. These issues can create run-to-run differences in the least significant digits. Two runs, same code, two answers. This is by design and is not an issue of correctness. Subtle changes in the memory layout of the program’s data, differences in loading of the ALU registers and operation order, and differences in threading all due to unrelated processes running on the same machine cause these run-to-run differences.

Managing Reproducibility

Most importantly, one should test code’s numerical results only to the precision that can be expected by the algorithm, input data, and finally the limits of floating point arithmetic. To do this in unit tests, compare floating point numbers carefully only to a fixed number of digits. The code snippet below compares two double numbers and returns true only if the numbers match to a specified number of digits.

private static bool EqualToNumDigits( double expected, double actual, int numDigits )
    {
      double max = System.Math.Abs( expected ) > System.Math.Abs( actual ) ? System.Math.Abs( expected ) : System.Math.Abs( actual );
      double diff = System.Math.Abs( expected - actual );
      double relDiff = max > 1.0 ? diff / max : diff;
      if ( relDiff <= DOUBLE_EPSILON )
      {
        return true;
      }
 
      int numDigitsAgree = (int) ( -System.Math.Floor( Math.Log10( relDiff ) ) - 1 );
      return numDigitsAgree >= numDigits;
    }

This type of comparison should be used throughout unit testing code. The full code listing, which we use for our internal testing, is provided at the end of this article.

If it is essential to enforce binary run-to-run reproducibility to the limits of precision, NMath provides a flag in its configuration class to ensure this is the case. However this flag should be set for unit testing only because there can be a significant cost to performance. In general, expect a 10% to 20% reduction in performance with some common operations degrading far more than that. For example, some matrix multiplications will take twice the time with this flag set.

Note that the number of threads that Intel’s MKL library uses ( which NMath depends on ) must also be fixed before setting the reproducibility flag.

int numThreads = 2;  // This must be fixed for reproducibility.
NMathConfiguration.SetMKLNumThreads( numThreads );
NMathConfiguration.Reproducibility = true;

This reproducibility run configuration for NMath cannot be unset at a later point in the program. Note that both setting the number of threads and the reproducibility flag may be set in the AppConfig or in environmental variables. See the NMath User Guide for instructions on how to do this.

Paul

References

M. A. Cornea-Hasegan, B. Norin. IA-64 Floating-Point Operations and the IEEE Standard for Binary Floating-Point Arithmetic. Intel Technology Journal, Q4, 1999.
http://gec.di.uminho.pt/discip/minf/ac0203/icca03/ia64fpbf1.pdf

D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic. Computing Surveys. March 1991.
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Full double Comparison Code

private static bool EqualToNumDigits( double expected, double actual, int numDigits )
    {
      bool xNaN = double.IsNaN( expected );
      bool yNaN = double.IsNaN( actual );
      if ( xNaN && yNaN )
      {
        return true;
      }
      if ( xNaN || yNaN )
      {
        return false;
      }
      if ( numDigits <= 0 )
      {
        throw new InvalidArgumentException( "numDigits is not positive in TestCase::EqualToNumDigits." );
      }
 
      double max = System.Math.Abs( expected ) > System.Math.Abs( actual ) ? System.Math.Abs( expected ) : System.Math.Abs( actual );
      double diff = System.Math.Abs( expected - actual );
      double relDiff = max > 1.0 ? diff / max : diff;
      if ( relDiff <= DOUBLE_EPSILON )
      {
        return true;
      }
 
      int numDigitsAgree = (int) ( -System.Math.Floor( Math.Log10( relDiff ) ) - 1 );
      //// Console.WriteLine( "x = {0}, y = {1}, rel diff = {2}, diff = {3}, num digits = {4}", x, y, relDiff, diff, numDigitsAgree );
      return numDigitsAgree >= numDigits;
    }

Special Functions

Monday, May 11th, 2015

Phi function, from Abramowitz & Stegun (1965) page 258. Motivated by the need of some special functions when writing signal processing code for NMath, we decided to add a suite of special functions to be included in NMath 6.1. While the field of special functions is vast [1], our 41 functions cover many of the most commonly needed functions in physics and engineering. This includes the gamma function and related functions, Bessel functions, elliptic integrals, and more. All special functions in NMath are now organized in a SpecialFunctions class, which is structured similarly to the existing StatsFunctions and NMathFunctions classes.

Special functions list

Below is a complete list of the special functions now available in the SpecialFunctions class which resides in the CenterSpace.NMath.Core name space. Previously a handful of these functions were available in either the NMathFunctions or StatsFunctions classes, but now those functions have been deprecated and consolidated into the SpecialFunctions class. Please update your code accordingly as these deprecated functions will be removed from NMath within two to three release cycles.

Using these special functions in your code is simple.

using namespace CenterSpace.NMath.Core
 
// Compute the Jacobi function Sn() with a complex argument.
var cmplx = new DoubleComplex( 0.1, 3.3 )
var sn = SpecialFunctions.Sn( cmplx, .3 );  // sn = 0.16134 - i 0.99834
 
// Compute the elliptic integral, K(m)
var ei = SpecialFunctions.EllipticK( 0.432 ); // ei = 1.80039

Below is a complete list of all NMath special functions.

 

Special Function Comments
EulerGamma A constant, also known as the Euler-Macheroni constant. Famously, rationality unknown.
Airy Provides solutions Ai, Bi, and derivatives Ai’, Bi’ to y” – yz = 0.
Zeta The Riemann zeta function.
PolyLogarithm The Polylogarithm, Li_n(x) reduces to the Riemann zeta for x = 1.
HarmonicNumber The harmonic number is a truncated sum of the harmonic series, closely related to the digamma function.
Factorial n!
FactorialLn The natural log of the factorial, ln( n! ).
Binomial The binomial coefficient, n choose k; The number of ways of picking k unordered outcomes from n possibilities.
BinomialLn The natural log of the binomial coefficient.
Gamma The gamma function, conceptually a generalization of the factorial.
GammaReciprocal The reciprocal of the gamma function.
IncompleteGammaFunction Computes the gamma integral from 0 to x.
IncompleteGammaComplement Computes the gamma integral from x to infinity (and beyond!).
Digamma Also known as the psi function.
GammaLn The natural log of the gamma function.
Beta The beta integral is also known as the Eulerian integral of the first kind.
IncompleteBeta Computes the beta integral from 0 to x in [0,1].
Ei The exponential integral.
Elliptic Integrals
EllipticK The complete elliptic integral, K(m), of the first kind. Note that m is related to the elliptic modulus k with, m = k * k.
EllipticE( m ) The complete elliptic integral, E(m), of the second kind.
EllipticF The incomplete elliptic integral of the first kind.
EllipticE(phi, m) The incomplete elliptic integral of the second kind.
EllipJ Computes the Jacobi elliptic functions Cn(), Sn(), and Dn() for real arguments.
Sn Computes the Jacobi elliptic function Sn() for complex arguments.
Cn Computes the Jacobi elliptic function Cn() for complex arguments.
Bessel Functions
BesselI0 Modified Bessel function of the first kind, order zero.
BesselI1 Modified Bessel function of the first kind, first order.
BesselIv Modified Bessel function of the first kind, non-integer order.
BesselJ0 Bessel function of the first kind, order zero.
BesselJ1 Bessel function of the first kind, first order.
BesselJn Bessel function of the first kind, arbitrary integer order.
BesselJv Bessel function of first kind, non-integer order.
BesselK0 Modified Bessel function of the second kind, order zero.
BesselK1 Modified Bessel function of the second kind, order one.
BesselKn Modified Bessel function of the second kind, arbitrary integer order.
BesselY0 Bessel function of the second kind, order zero.
BesselY1 Bessel function of the second kind, order one.
BesselYn Bessel function of the second kind of integer order.
BesselYv Bessel function of the second kind, non-integer order.
Hypergeometric Functions
Hypergeometric1F1 The confluent hypergeometric series of the first kind.
Hypergeometric2F1 The Gauss or generalized hypergeometric function.

 

Let us know if you need any additional special functions and we’ll see if we can add them.

Mathematically,

Paul Shirkey

References

[1] Abramowitz, M. and Stegun, I. (1965). Handbook of Mathematical Functions. Dover Publications. ( Abramowitz and Stegun PDF )
[2] Wolfram Alpha LLC. (2014). www.wolframalpha.com
[3] Weisstein, Eric W. “[Various Articles]” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/
[4] Moshier L. Stephen. (1995) The Cephes Math Library. (Cephes).

Announcing NMath 6.1 and NMath Stats 4.1

Thursday, April 30th, 2015

We’re pleased to announce new versions of the NMath libraries – NMath 6.1, and NMath Stats 4.1.

Added functionality includes:

  • Upgraded to Intel MKL 11.2 Update 2 with resulting performance increases.
  • Updated NMath Premium GPU code to CUDA 6.
  • Added classes for solving linear and nonlinear programming problems with integer or binary constraints.
  • Added class SpecialFunctions containing special functions such as factorial, binomial, the gamma function and related functions, Bessel functions, elliptic integrals, and many more. (Prior versions of a few of these functions, such as StatsFunctions.IncompleteGamma, are now deprecated.)
  • Added a new native library, nmath_sf_x86.dll and nmath_sf_x64.dll, with high-performance C language implementations of the special functions.
  • Added single-precision versions of general sparse matrix and vector types.

For more complete changelogs, see here and here.

Upgrades are provided free of charge to customers with current annual maintenance contracts. To request an upgrade, please visit our upgrade page, or contact sales@centerspace.net. Maintenance contracts are available through our webstore.

Announcing NMath 6.0 and NMath Stats 4.0

Tuesday, August 19th, 2014

We’re pleased to announce new versions of the NMath libraries – NMath 6.0, and NMath Stats 4.0.

Added functionality includes:

  • Upgraded to Intel MKL 11.1 Update 3 with resulting performance increases.
  • Added Adaptive Bridge™ technology to NMath Premium edition, with support for multiple GPUs, per-thread control for binding threads to GPUs, and automatic performance tuning of individual CPU–GPU routing to insure optimal hardware usage.
  • NMath linear programming, nonlinear programming, and quadratic programming classes are now built on the Microsoft Solver Foundation (MSF). The Standard Edition of MSF is included with NMath.
  • Added classes for solving nonlinear programming problems using the Stochastic Hill Climbing algorithm, for solving quadratic programming problems using an interior point algorithm, and for solving constrained least squares problems using quadratic programming methods.
  • Added support for MKL Conditional Numerical Reproducibility (CNR).

For more complete changelogs, see here and here.

Upgrades are provided free of charge to customers with current annual maintenance contracts. To request an upgrade, please contact sales@centerspace.net. Maintenance contracts are available through our webstore.