Blog

Archive for the ‘NMath’ Category

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.

Optimal Portfolio Allocation

Wednesday, January 22nd, 2014

The problem of optimal portfolio allocation, in its simplest form, asks the question of how to fully allocate a given amount of wealth across a fixed universe of investments to achieve a minimum-risk goal-expected return. The known quantities are the potential field of investments, their performance history, and the goal rate of return; The unknown is the wealth allocation across the investments. In this standard formulation, the phrase ‘minimum-risk’ is understood to be equivalent to ‘minimum-variance’ of the portfolio return. Said another way, in this model the investor will expect to get a ‘risk-premium’ (higher return) for investing in assets which have a higher historic return variance. This understanding of risk is at the heart of Harry Markowitz’s portfolio theory [5].

In the following brief development of optimal portfolio allocation theory we follow the paper by P. J. Atzberger [2]. If we have n investments we represent a fully invested portfolio as,


\sum^{n}_{i=1}w_i = 1

where the fraction of wealth invested in the ith asset is wi. We use this normalization to avoid working with the absolute magnitudes of the assets and resultant portfolios.

The rate of return of the entire portfolio, ρp at a given time t is,


\rho_p(t) = \sum^{n}_{i=1}w_i*\rho_i

or simply the weighted sum of investment returns.

And because expectation is a linear operator, the expected rate of return of a given portfolio is,


\mu_p = E \left( \sum^{n}_{i=1}w_i*\rho_i \right ) = \sum^{n}_{i=1}w_i*\mu_i

Where μi is expected rate of return of the ith investment option.

And finally from the expected return of return, the portfolio variance will be given by,


\sigma_p^2 = E \left( |\rho_p - \mu_p|^2 \right ) = \textbf{w}^T V \textbf{w}

where w is the vector of normalized fractional investments defined above. (See eq. (13) in [2] for the full derivation). This is the minimum theoretical framework needed to setup the optimization problem.

Optimization with Markowitz Theory

Using our established problem framework we can cast the optimal portfolio allocation problem as a constrained quadratic optimization problem by minimizing the portfoilo’s variance (risk) subject to two linear constraints.


min {\,\,\,} {1\over{2}} \sum^{n}_{i,j=1} \omega_i \omega_j \sigma_{i,j}

subject to,


\sum^{n}_{i=1} \omega_i \mu_j - \mu_p = 0

and,


\sum^{n}_{i=1} \omega_i - 1 = 0

The minimized objective function is the variance of the portfolio – which is thought of as risk equivalent in this model. The first constraint insures that the built portfolio achieves a goal return of μp and the second constraint enforces the full allocation of wealth. There are many variants of this basic setup such as including a zero-risk investment (cash) to the pool of possible investments or adding constraints to allocate various percentages of wealth into different risk-classified investments [2][4].

Solving the Quadratic Programming Problem

This constrained optimization problem can be solved using NMath‘s quadratic programming class, ActiveSetQPSolver(). The basic steps of setting up the QP solver involve (1) Reading and parsing the investment data, (2) Setting up the problem constraints, (3) Executing the solve. The data for this example come from a data repository for exploring various Operations Research (OR) problems, run by J.E. Beasley. The specific data we study here has a universe of 225 potential assets, is found here, and is the largest test data set used in the paper “Heuristics for cardinality constrained portfolio optimisation” [3] (A brief description of the data set can be found on the same server here). The data set provides the basic information we need for optimizing the portfolio – the mean return and standard deviation of return of each investment and the correlation between all possible pairs of assets. The code to parse this file is appended to the end of this article.

Before generating the QP’s constraints we need to set up a QuadraticProgrammingProblem instance.

  // Set up QR problem using port5.txt data
  FloatMatrix covariance = correlation * 
                           ( NMathFunctions.Product( new FloatMatrix( stdDevReturn ), 
                           ( new FloatMatrix( stdDevReturn ).Transpose() ) ) );
  QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem();
  problem.H = covariance;
  problem.c = new DoubleVector(n, 0);

First, the covariance matrix is computed from the investment standard deviations and correlations. Next we create and initialize a QuadraticProgrammingProblem instance. Since this QP solver class solves programs of the form: Minimize x'Hx + x'c, the c variable is set to zero.

The following code snippet sets up the QP constraints.

  ...
  QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem();
  ...
  // Set up constraints
  LinearConstraint sumToOne = new LinearConstraint( 
                   new DoubleVector( n, 1 ), 1, ConstraintType.EqualTo );
  LinearConstraint returnRate = new LinearConstraint( 
                   meanReturn, targetReturn, ConstraintType.GreaterThanOrEqualTo );
  problem.AddConstraint( returnRate );
  problem.AddConstraint( sumToOne );
 
  // Add the bound on each var, w: w > 0
  for ( int i = 0; i < n; i++ )
    problem.AddLowerBound( i, 0.0 );

The first two constraints, sumToOne and returnRate were outlined above in the brief Markowitz theory section. The third constraint has been added to insure that we have a positive ownership in all assets; in this way we have disallowed shorting assets. This constraint can be deleted to include shorted asset in the optimized portfolio, if desired. All of these constraints are added to the QuadraticProgrammingProblem instance.

We can now put all of this together and solve the constrained 225 variable QP using the ActiveSetQPSolver class [0].

Asset   Allocation
   8,   0.0795
  39,   0.0866
  42,   0.0812
  59,   0.1201
  61,   0.2567
  96,   0.0593
 128,   0.0741
 170,   0.0573
 195,   0.0980
 214,   0.0688
 224,   0.0183
(All other assets have an allocation of zero.)
 
Total allocation:  1.0000
Total return-rate:  0.2000%, Target return-rate  0.2000%
Total solver time required:  0.1830 seconds

The optimal allocation includes 11 of our 225 investments and closely adheres to all constraints. The complete code listing is provided at the end of this article.

Happy Computing,

Paul Shirkey

References, Resources, Notes

[0] Due to the size of this QP problem and the number of constraints (227) the current release of NMath‘s ActiveSetQPSolver() fails to correctly find the optimal portfolio. However, the next release of NMath due in May (NMath 6.0) includes a major upgrade to many of our LP, QP, and NLP classes. Many of these classes in NMath 6.0 are backed by the Microsoft Solver Foundation, providing the ability to solve much larger problems correctly and efficiently.

[1] The port5 data set can be found at http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/port5.txt . This data set is part of the OR-Library, a set of standard data sets maintained by J. E. Beasley at http://people.brunel.ac.uk/~mastjjb/jeb/info.html.

[2] This blog note follows the notation of this clearly written paper, by Paul J. Atzberger, on optimal portfolio allocation, variants, and the efficient frontier: http://www.math.ucsb.edu/~atzberg/finance/portfolioTheory.pdf.

[3] Chang, T.-J., Meade, N., Beasley, J.E. and Sharaiha, Y.M., “Heuristics for cardinality constrained portfolio optimisation” Comp. & Opns. Res. 27 (2000) 1271-1302.

[4] MATLAB uses this same data set in their optimization documentation.

[5] Markowitz, H.M. (March 1952). “Portfolio Selection”. The Journal of Finance 7 (1): 77–91. [paper in pdf]

.NET Code Listings

Full QP solution code. Note that NMath 6.0 (due May 2014) is required.

using System;
using CenterSpace.NMath.Core;  
using CenterSpace.NMath.Analysis;
using CenterSpace.NMath.Charting.Microsoft;
using System.Windows.Forms.DataVisualization.Charting;
 
public void OptimalPortAllocation()
    {
      // Get port5 data set : http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/port5.txt
      // File Format of port5:
      // Line 1: Text description of data
      // Line 2: number of assets, N
      // Line 2 - 2+N : list of mean return _ standard return
      // Remaining 1/2 N*N lines: correlation matrix between asset i and j
      FloatVector meanReturn = null ;
      FloatVector stdDevReturn = null;
      FloatMatrix correlation = null;
      int n = loadData( ref meanReturn, ref stdDevReturn, ref correlation );
 
      double targetReturn = 0.002;
 
      // Set up QR problem using port5.txt data
      FloatMatrix covariance = correlation * ( NMathFunctions.Product( new FloatMatrix( stdDevReturn ), ( new FloatMatrix( stdDevReturn ).Transpose() ) ) );
      QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem();
      problem.H = covariance;
      problem.c = new DoubleVector(n, 0);
 
      // Set up constraints
      LinearConstraint sumToOne = new LinearConstraint( new DoubleVector( n, 1 ), 1, ConstraintType.EqualTo );
      LinearConstraint returnRate = new LinearConstraint( meanReturn, targetReturn, ConstraintType.GreaterThanOrEqualTo );
      problem.AddConstraint( returnRate );
      problem.AddConstraint( sumToOne );
 
      // Add the bound on each var, w: w > 0
      for ( int i = 0; i < n; i++ )
        problem.AddLowerBound( i, 0.0 );
 
      // Solve constrainted QP
      var solver = new ActiveSetQPSolver();
      solver.StepSizeEpsilon = 1e-10;
 
      var startEqualAllocation = new DoubleVector( n, 1.0 / n );
      var optAllocation = new DoubleVector();
 
      long before = DateTime.Now.Ticks;
      if ( solver.Solve( problem, startEqualAllocation ) )
      {
        optAllocation = solver.OptimalX;
      }
      else
      {
        throw new Exception( "Solve failed to find a solution." );
      }
      long after = DateTime.Now.Ticks;
      long ticks = after - before;
      double millis = ( (double) ticks ) / 10000;
 
      Console.WriteLine("Asset #    Allocation");
      double wealth = 0.0;
      double measuredReturnRate = 0.0;
      for ( int i = 0; i < optAllocation.Length; i++ )
      {
        if ( optAllocation[i] > 0.0000001 )
        {
          Console.WriteLine( String.Format( "{0,4},  {1,7:0.0000}",i, optAllocation[i]));
          wealth += optAllocation[i];
          measuredReturnRate += optAllocation[i] * meanReturn[i];
        }
 
      }
      Console.WriteLine(String.Format("Total allocation: {0,7:0.0000}", wealth) );
      Console.WriteLine( String.Format( "Total return rate: {0,7:0.0000}%, Target return rate {1,7:0.0000}%", 100 * measuredReturnRate, 100 * targetReturn ) );
      Console.WriteLine( String.Format( "Total solver time required: {0,7:0.0000} seconds", millis/1000 ) );
 
 
      // Chart optimal portfolio allocation
      DoubleVector x = new DoubleVector(n, 0, 1);
      DoubleVector w = optAllocation;
      var chart = NMathChart.ToChart( x, w );
      chart.Series[0].MarkerStyle = MarkerStyle.Square;
      chart.Series[0].ChartType = SeriesChartType.Bar;
      NMathChart.Show( chart );
 
    }

Below is the full parsing code for the port5.txt data file.

    private int loadData( ref FloatVector meanReturn, ref FloatVector stdDevReturn, ref FloatMatrix correlation )
    {
      int n;
 
      // Fix up the path here to point to your data file copy.
      using ( StreamReader file = new StreamReader( "./data/PortfolioOptimzation_port5.data.txt" ) )
      {
        file.ReadLine(); // skip first line of explanatory text
        n = Int32.Parse( file.ReadLine() );
 
        // Now get the next n mean return and standard return values.
        meanReturn = new FloatVector( n );
        stdDevReturn = new FloatVector( n );
        string[] values = new string[3];
        for ( int k = 0; k < n; k++ )
        {
          values = file.ReadLine().Split( ' ' );
          meanReturn[k] = float.Parse( values[1] );
          stdDevReturn[k] = float.Parse( values[2] );
        }
 
        // Now get the correlation matrix.
        correlation = new FloatMatrix( n, n );
        int cnt = 0;
        for ( int i = cnt; i < n; i++ )
        {
          for ( int j = cnt; j < n; j++ )
          {
            values = file.ReadLine().Split( ' ' );
            correlation[i, j] = float.Parse( values[3] );
            correlation[j, i] = float.Parse( values[3] );
          }
 
          cnt++;
        }
      }
      return n;
    }

Chebyshev Filters in C#

Wednesday, December 11th, 2013

There are three classes of widely used IIR (recursive) filters in signal processing: Butterworth, Chebyshev, and elliptical. In this article I will give a short introduction to the Chebyshev filter, present a code implementation, and end with a usage example. The Butterworth filter was discussed in a previous blog article.

Chebyshev filters come in two flavors defined by either allowing ripple in the pass-band (type 1) or ripple the stop-band (type 2). We will only discuss the type 1 filters here, as the type 2 filters are rarely used. Type 1 Chebyshev filters trade-off steeper roll-off with ripple in the pass band and in the limit of zero ripple, mimic Butterworth filters. Simply put, the more ripple allowed in the passband, the steep roll-off the filter achieves.
Chebyshev filters derive their name from Chebyshev polynomials [1]. Finding Chebyshev polynomials for a given order represent a computational challenge, and when computed often exhibit various numerical instabilities (often overflow) if not done properly. For an excellent review and comparison of computational techniques for computing Chebyshev polynomials see this paper by Wolfram Koepf. Here, we will side step the entire computational issue, as legions of engineers have done before us, and simply look up the Chebyshev polynomials, by order, in a switch statement. Because Chebyshev filters with an order n above 10 or 12 would be unusual, this table look-up approach can be considered, practically speaking, complete.

Wikipedia has a excellent overview article discussing the Chebyshev filter, so I’ll refer you there for a background treatment of Chebyshev filtering theory.

Chebyshev Filters

The gain expression for Chebyshev filters has a very similar structure to Butterworth filters. As with the Butterworth filter, there are the three familiar design parameters, DC Gain, cutoff frequency f0, and filter order n, with an additional parameter, epsilon, for specifying the desired pass-band ripple.


Chebyshev gain equation

Chebyshev gain equation


The Tn signifies a Chebyshev polynomial of order n. So for example when n = 2, it happens that Tn = 2*x*x - 1, where when substituted into the gain equation we replace x with f/f0. The ripple factor can be equated to pass-band ripple in Db with the equation:


Passband ripple equation expressed in decibles

Passband ripple in Db

Typical values for epsilon would range between 0.1 and 1.

The Chebyshev Filter in Code

We take the identical approach to implementing the Chebyshev filter in code as we did with the Butterworth filter. We will first compute the input signal’s FFT, then multiply that by the above filter gain, and then take the inverse FFT of that product resulting in our filtered signal. This is a O( n*log(n)) operation.

But before we do that, we need a way to find the Chebyshev polynomials of a given order n. To that end I’ve build a simple method that returns a Chebyshev polynomial given the order; if the requested order is larger than the look-up table, a NotImplementedException is thrown. This look-up approach can easily be extended as far as would be needed for most any practicable Chebyshev filter.

Chebyshev Polynomial Lookup

 private Polynomial ChebyshevPolynomial( int order)
    {
      switch( order )
      {
        case 0:
          return new Polynomial( new DoubleVector( "1" ) );
        case 1:
          return new Polynomial( new DoubleVector( "0 1" ) );
        case 2:
          return new Polynomial( new DoubleVector( "-1 0 2" ) );
        case 3:
          return new Polynomial( new DoubleVector( "0 -3 0 4" ) );
        case 4:
          return new Polynomial( new DoubleVector( "1 0 -8 0 8" ) );
        case 5:
          return new Polynomial( new DoubleVector( "0 5 0 -20 0 16" ) );
        case 6:
          return new Polynomial( new DoubleVector( "-1 0 18 0 -48 0 32" ) );
        case 7:
          return new Polynomial( new DoubleVector( "0 7 0 56 0 -112 0 64" ) );
        case 8:
          return new Polynomial( new DoubleVector( "1 0 -32 0 160 0 -256 0 128" ) );
        default:
          throw new NotImplementedException();     
      }
   }

An excellent on-line resource for Chebyshev polynomials, or any other integer sequence can be found at OEIS, the Online Encyclopaedia of Integer Sequences. If you wish to delve into computing the coefficients of Chebyshev polynomials, a good review of techniques can be found here.

Chebyshev Filter

The following code implements the Chebyshev filter, following the three steps mentioned above: [1] FFT(signal), [2] FFT(signal)*G(f), [3] IFFT.

private DoubleComplexVector ChebeshevFilter( DoubleComplexVector signal, 
  double sampleFrequency, int order, double f0, 
  double ripple, double DCGain )
    {
      // Get the Chebyshev polynomial
      Polynomial T = ChebyshevPolynomial( order );
 
      int N = signal.Length;
 
      // Apply forward FFT in place
      DoubleComplexForward1DFFT fft = new DoubleComplexForward1DFFT( N );
      DoubleComplexVector signalFFT = fft.FFT( signal );
 
      // Remove DC offset
      signalFFT[0] = 0;
 
      if ( f0 > 0 )
      {
 
        double numBins = N / 2;  // Half the length of the FFT by symmetry
        double binWidth = sampleFrequency / numBins; // Hz
 
        // Filter
        System.Threading.Tasks.Parallel.For( 1, N / 2, i =>
        {
          double binFreq = binWidth * i;
          double gain = 
            DCGain / ( Math.Sqrt( ( 1 + Math.Pow( T.Evaluate( binFreq / f0) , 2.0 ) ) ) );
          signalFFT[i] *= gain;
          signalFFT[N - i] *= gain;
        } );
 
      }
 
      // Reverse filtered signal
      DoubleComplexBackward1DFFT ifft = new DoubleComplexBackward1DFFT( N );
      ifft.SetScaleFactorByLength(); // Needed to get the correct amplitude
      signal = ifft.FFT( signalFFT );
 
      return signal;
    }

Filtering Example

To test the Chebyshev filtering code above, I’ve generated a test signal which is the sum to two sinusoids of frequences 1 Hz and 8 Hz. For comparison, this is the same test signal I used to test the Butterworth filter in my previous post. The test signal is sampled at 25 Hz and then filtered using a third order Chebyshev filter with a frequency cut-off of 1.5 Hz and ripple factor (epsilon) of 0.2.

public void ChebyshevFilterTest()
    {
      DoubleComplexVector signal = new DoubleComplexVector( 100 );
 
      double Fs = 25; // Sample rate, 25Hz
      double f1 = 1; // 1 Hz Signal
      double f2 = 8; // 8 Hz Signal
      for ( int t = 0; t < 100; t++ )
      {
        signal[t] = new DoubleComplex( Math.Sin( 2 * Math.PI * f1 / Fs * t ) + Math.Sin( 2 * Math.PI * f2 / Fs * t ), 0.0 );
      }
 
      // Filter      
      DoubleComplexVector filteredSignal = this.ChebeshevFilter( signal, 25, 3, 1.5, 0.2, 1.0 );
 
      // Display results
      Chart chart = NMathChart.ToChart( new DoubleVector( signal.ToRealDataTable() ) );
      chart.Titles[0].Text = "Input signal, sampled at 25 Hz";
      chart.ChartAreas[0].AxisY.Title = "Input signal magnitude";
      chart.ChartAreas[0].AxisX.Minimum = 0;
      Chart chartfft = NMathChart.ToChart( new DoubleVector( filteredSignal.ToRealDataTable() ) );
      chartfft.Titles[0].Text = "Filtered signal";
      chartfft.ChartAreas[0].AxisY.Title = "Filtered signal magnitude";
      chartfft.ChartAreas[0].AxisX.Minimum = 0;
 
      List<Chart> charts = new List<Chart>() { chart, chartfft };
 
      Chart compositeChart = NMathChart.Compose( charts, 1, 2, NMathChart.AreaLayoutOrder.ColumnMajor );
      compositeChart.Width = 650;
      compositeChart.Height = 300;
 
      NMathChart.Show( compositeChart );
 
    }

This example test code, using NMath charts, generates the following chart of the unfiltered and filtered signal side by side.

Chebyshev filtering example

-Happy Computing,

Paul

[1] Chebyshev is a transliteration of a Russian name and in the past has been variously spelled as Chebychev, Tschebyscheff, Tchebysheff and Tchebichef.

NMath Tutorial videos

Thursday, October 31st, 2013

We are proud to announce a series of NMath tutorial videos on how to use CenterSpace’s math library products. We are starting, naturally, with Getting Started with NMath.

You can download it here: MP4

Please let us know which topics you want us to cover. Email support@centerspace.net

Cheers,
Trevor

Butterworth Filters in C#

Wednesday, October 30th, 2013

There are three classes of widely used IIR (recursive) filters in signal processing: Butterworth, Chebyshev, and elliptical. In this article I will discuss the Butterworth filter and provide example code implementing and using the filter. The Chebyshev and elliptical filters will be discussed in follow up articles.

Butterworth filters are desirable for their ease of implementation, good phase response, and their smooth monotonic frequency response in both the pass-band and the stop-band. However they don’t achieve the steep roll-off of Chebyshev filters, but Chebyshev filters achieve this steeper roll-off by allowing either some ripple in the pass-band (type 1) or the stop-band (type 2). The type 2 Chebyshev filters are rarely used so only type 1 Chebyshev filters will be discussed in follow up articles. The elliptical filters strike a balance between the type 1 and type 2 Chebyshev filters by allowing ripple in both the stop and pass-bands while achieving the steepest roll-off gain of all three filter types. Wikipedia has excellent overview articles discussing the Butterworth, the Chebyshev, and the elliptical filters, so I’ll refer you to those articles for a thorough background treatment of each filter type instead of repeating that information here.

Butterworth Filters

Butterworth filters exhibited a ripple free frequency response with a -20*n Db/decade roll-off at the cutoff frequency, where n is the order of the filter. There are only three design parameters for a Butterworth filter, the order n, the cut-off frequency Cutoff frequency, and the DC gain, DC gain, or the gain at zero frequency. So the gain of any Butterworth filter can be written in terms of these three parameters.


Butterworth filter gain

Butterworth filter gain


Engineers typically think in terms of Hertz, so we can relate the linear frequency cutoff we are accustom to, to the angular frequency we have above, with the familiar equation, Angular cut-off frequency. Now rewriting the above equation in terms of frequency we have:


Butterworth filter gain in terms of f

Butterworth filter gain in terms of f

The Butterworth Filter in C#

To digitally implement a Butterworth filter in code, we work in the frequency domain (for performance reasons), following the three steps outlined in the pseudo code below.

function ButterworthFilter
{
  fftSignal = FFT( signal );
 
  // Scale FFT by butterworth gain
  for( int w = 0; w < N/2; w++)
    filteredSignal[w] = fftSignal[w] * butterworthGain[w];
 
  return inverseFFT( filteredSignal );
}

Given the dual relationship between convolution and multiplication in the time and frequency domains, we could equally well implement this filter by convolving in the time domain the inverse fourier of the Butterworth’s frequency response with the signal, but this would require O(n*n) time. By using the FFT, and working in the frequency domain, we achieve the same effect in O(n* log(n)) time, where n is the length of the signal to be filtered.

DoubleComplexVector ButterworthFilter( DoubleComplexVector signal, 
   double sampleFrequency, int order, double f0, double DCGain )
    {
 
      var N = signal.Length;
 
      // Apply forward FFT
      var fft = new DoubleComplexForward1DFFT( N );
      var signalFFT = fft.FFT( signal );
 
      if ( f0 > 0 )
      {
 
        var numBins = N / 2;  // Half the length of the FFT by symmetry
        var binWidth = sampleFrequency / numBins; // Hz
 
        // Filter
        System.Threading.Tasks.Parallel.For( 1, N / 2, i =>
        {
          var binFreq = binWidth * i;
          var gain = DCGain / ( Math.Sqrt( ( 1 + 
                        Math.Pow( binFreq / f0, 2.0 * order ) ) ) );
          signalFFT[i] *= gain;
          signalFFT[N - i] *= gain;
        } );
 
      }
 
      // Reverse filtered signal
      var ifft = new DoubleComplexBackward1DFFT( N );
      ifft.SetScaleFactorByLength(); // Needed to get the correct amplitude
      signal = ifft.FFT( signalFFT );
 
      return signal;
   }

This general Butterworth filtering code begs a bit of explanation. The FFT provides the Fourier transform for N/2 positive frequencies, and symmetrically (via a complex conjugation) for N/2 negative frequencies. Because of this we can loop just N/2 times scaling the symmetric halves simultaneously by our Butterworth gain in our Parallel.For loop for improved performance. When scaling the FFT by the gain of the Butterworth filter, we must know the angular frequency of each bin we scale, which is captured by the binFreq variable; recalling that the frequency resolution of a length N FFT is ( Fs / N / 2 ) Hz, where Fs is the sampling frequency.

Although this code will execute very efficiently it can be significantly optimized if it is specialized to filtering only signals of a fixed length and fixed sampling frequency – a very common situation in many systems. In this case the Butterworth gains can all be computed ahead, avoiding the expensive Math.Sqrt() calls, also the NMath FFT classes would only need to be created once outside the filtering routine. Finally we could do all of the computation in-place and eliminate the copy to a new result vector, with:

   ...
   // Apply forward FFT in place
   var fft = new DoubleComplexForward1DFFT( N );
   fft.FFTInPlace( signal );
   ...

Filtering Example

To test our new Butterworth filter, I generated a test signal which is simply the sum of two sinusoids of frequencies 1 Hz and 8 Hz. This continuous signal is then sampled at 25 Hz and filtered with a third order Butterworth filter using a cut-off frequency of 1.5 Hz and a DC gain of one. Finally, a chart is generated of the input signal and of the filtered signal to visualize the action of the filter in the time domain.

    public void FilterTest()
    {
      var signal = new DoubleComplexVector( 100 );
 
      var Fs = 25; // Sample rate, 25Hz
      var f1 = 1; // 1 Hz Signal
      var f2 = 8; // 8 Hz Signal
      for ( int t = 0; t < 100; t++ )
      {
        signal[t] = new DoubleComplex( Math.Sin( 2 * Math.PI * f1 / Fs * t ) + 
                    Math.Sin( 2 * Math.PI * f2 / Fs * t ), 0.0 );
      }
 
      // Filter      
      var filteredSignal = ButterworthFilter( signal, 25, 3, 1.5, 1.0 );
 
      // Display results
      var chart = NMathChart.ToChart( new DoubleVector( signal.ToRealDataTable() ) );
      chart.Titles[0].Text = "Input signal, sampled at 25 Hz";
      chart.ChartAreas[0].AxisY.Title = "Input signal magnitude";
      var chartfft = NMathChart.ToChart( new DoubleVector( filteredSignal.ToRealDataTable() ) );
      chartfft.Titles[0].Text = "Filtered signal";
      chartfft.ChartAreas[0].AxisY.Title = "Filtered signal magnitude";
 
      List<Chart> charts = new List<Chart>() { chart, chartfft };
 
      var compositeChart = NMathChart.Compose( charts, 1, 2, NMathChart.AreaLayoutOrder.ColumnMajor );
      compositeChart.Width = 650;
      compositeChart.Height = 300;
      NMathChart.Show( compositeChart );
 
    }

This example code generates the following composite chart showing on the left the unfiltered signal, and on the right the filtered signal.

Butterworth filtering example

Butterworth filtering example

As expected, almost all of the high frequency sinusoid at 8Hz has been eliminated from the original signal. In this case where the ‘noisy’ 8 Hz signal is well separated from our 1Hz signal, the Butterworth filter does an excellent job. The next signal processing blog entry will cover the Chebyshev type 1 filter.

-Happy Computing,

Paul Shirkey

Please contact us at sales@centerspace.net with any questions you may have.