Blog

Archive for the ‘NMath’ Category

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.

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.