Adds Remez filter designer for low, high and band-pass filters. Support

for Hilbert and Differentiator FIR filter designers may be added in the
future.

Adds developer filter viewer tool written in JavaFX to visualize the
frequency and impulse response of FIR filters.  This application is not
intended to be a user application and is specifically for the developer
to analyze filters and filter specifications that will be used in the
code.

NOTE: Java 8 JDK is now required to build this project for version 0.3.x
and forward.  The JavaFX visualizer application and future user
interfaces will use the Java 8 version of JavaFX

.classpath

@@ -1,7 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <classpath>
 	<classpathentry kind="src" path="src"/>
-	<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER/org.eclipse.jdt.internal.debug.ui.launcher.StandardVMType/JavaSE-1.7"/>
 	<classpathentry kind="lib" path="imports/commons-lang3-3.2.1.jar"/>
 	<classpathentry kind="lib" path="imports/jide-oss-3.5.15.jar"/>
 	<classpathentry kind="lib" path="imports/jsr80-1.0.1.jar"/>
@@ -30,5 +29,6 @@
 			<attribute name="javadoc_location" value="jar:file:/home/denny/Downloads/tablefilter-swing-5.2.1-javadoc.jar!/"/>
 		</attributes>
 	</classpathentry>
+	<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER"/>
 	<classpathentry kind="output" path="bin"/>
 </classpath>

.project

@@ -5,6 +5,11 @@
 	<projects>
 	</projects>
 	<buildSpec>
+		<buildCommand>
+			<name>org.eclipse.xtext.ui.shared.xtextBuilder</name>
+			<arguments>
+			</arguments>
+		</buildCommand>
 		<buildCommand>
 			<name>org.eclipse.jdt.core.javabuilder</name>
 			<arguments>
@@ -13,5 +18,6 @@
 	</buildSpec>
 	<natures>
 		<nature>org.eclipse.jdt.core.javanature</nature>
+		<nature>org.eclipse.xtext.ui.shared.xtextNature</nature>
 	</natures>
 </projectDescription>

src/dsp/filter/design/FilterDesignException.java

@@ -0,0 +1,16 @@
+package dsp.filter.design;
+
+public class FilterDesignException extends Exception
+{
+    private static final long serialVersionUID = 1L;
+
+    public FilterDesignException( String text )
+    {
+        super( text );
+    }
+    
+    public FilterDesignException( String text, Throwable throwable )
+    {
+        super( text, throwable );
+    }
+}

src/dsp/filter/design/FilterView.css

@@ -0,0 +1,39 @@
+/* ....Show License.... */
+.chart {
+    -fx-background-image: url("../../../../images/Curve-fitted-background.png");
+    -fx-padding: 15 25 15 15;
+}
+.chart-plot-background {
+    -fx-background-color: red;
+    -fx-background-image: url("../../../../images/Curve-fitted-chart-background.png"), 
+                url("../../../../images/Curve-fitted-graph-gridlines.png");
+    -fx-background-size: cover, auto;
+    -fx-background-repeat: no-repeat, repeat;
+    -fx-background-position: 0% 0%, 0% 100%;
+    -fx-border-color: black black transparent transparent;
+}
+.chart-area-symbol {
+    -fx-background-color: white;
+}
+.chart-series-area-line { 
+    -fx-stroke: white; 
+    -fx-stroke-width: 2px; 
+}
+.chart-series-area-fill { 
+    -fx-fill: linear-gradient(to right, white, rgba(255,255,255,0)); 
+    -fx-blend-mode: OVERLAY;
+}
+.axis {
+    -fx-tick-mark-visible: false;
+    -fx-minor-tick-visible: false;
+    -fx-tick-length: 3;
+    -fx-minor-tick-length: 0;
+    -fx-border-color: transparent;
+    -fx-tick-label-fill: white;
+}
+.axis:bottom {
+    -fx-border-color: black transparent transparent transparent;
+}
+.axis:left {
+    -fx-border-color: transparent black transparent transparent;
+}

src/dsp/filter/design/FilterView.java

@@ -0,0 +1,171 @@
+package dsp.filter.design;
+
+import org.jtransforms.fft.FloatFFT_1D;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import dsp.filter.FilterFactory;
+import dsp.filter.Window.WindowType;
+import dsp.filter.fir.FIRFilterSpecification;
+import dsp.filter.fir.remez.RemezFIRFilterDesigner;
+import gui.control.CurveFittedAreaChart;
+import javafx.collections.ObservableList;
+import javafx.scene.chart.NumberAxis;
+import javafx.scene.chart.XYChart;
+import javafx.scene.chart.XYChart.Series;
+import javafx.scene.layout.BorderPane;
+import javafx.scene.layout.HBox;
+import javafx.scene.layout.Priority;
+
+public class FilterView extends BorderPane
+{
+	private final static Logger mLog = LoggerFactory.getLogger( FilterView.class );
+
+	private CurveFittedAreaChart mImpulseChart;
+	private NumberAxis mImpulseXAxis;
+	private NumberAxis mImpulseYAxis;
+	private CurveFittedAreaChart mDFTChart;
+	private NumberAxis mDFTXAxis;
+	private NumberAxis mDFTYAxis;
+	private static final int FFT_SIZE = 4096;
+	private FloatFFT_1D mFFT = new FloatFFT_1D( FFT_SIZE );
+	
+	public FilterView( float[] taps )
+	{
+		init();
+
+		if( taps != null )
+		{
+			setImpulseTaps( taps );
+
+			float[] dft = new float[ FFT_SIZE ];
+			System.arraycopy( taps, 0, dft, 0, taps.length );
+			
+			mFFT.realForward( dft );
+			
+			setDFTTaps( convertDFTToDecibel( dft, taps.length ) );
+		}
+	}
+	
+	private void init()
+	{
+		HBox hbox = new HBox();
+		
+		mImpulseXAxis = new NumberAxis( 0.0, 0.5, 1.0 );
+		mImpulseYAxis = new NumberAxis( -0.2, 1.2, .1 );
+		
+		mImpulseChart = new CurveFittedAreaChart( mImpulseXAxis, mImpulseYAxis );
+		mImpulseChart.setLegendVisible( false );
+		mImpulseChart.setHorizontalGridLinesVisible( false );
+		mImpulseChart.setVerticalGridLinesVisible(false);
+        mImpulseChart.setAlternativeColumnFillVisible(false);
+        mImpulseChart.setAlternativeRowFillVisible(false);
+        
+        String css = FilterView.class.getResource( "FilterView.css" ).toExternalForm();
+        mImpulseChart.getStylesheets().add( css );
+
+        hbox.getChildren().add( mImpulseChart );
+        HBox.setHgrow( mImpulseChart, Priority.ALWAYS );
+        
+		mDFTXAxis = new NumberAxis( 0.0, FFT_SIZE / 2, FFT_SIZE / 8 );
+		mDFTYAxis = new NumberAxis( -200.0, 10.0, 10.0 );
+		
+		mDFTChart = new CurveFittedAreaChart( mDFTXAxis, mDFTYAxis );
+		mDFTChart.setLegendVisible( false );
+		mDFTChart.setHorizontalGridLinesVisible( false );
+		mDFTChart.setVerticalGridLinesVisible(false);
+        mDFTChart.setAlternativeColumnFillVisible(false);
+        mDFTChart.setAlternativeRowFillVisible(false);
+        mDFTChart.getStylesheets().add( css );
+        
+        hbox.getChildren().add( mDFTChart );
+        HBox.setHgrow( mDFTChart, Priority.ALWAYS );
+
+        setCenter( hbox );
+	}
+	
+	public void setImpulseTaps( double[] taps )
+	{
+		float[] convertedTaps = new float[ taps.length ];
+		
+		for( int x = 0; x < taps.length; x++ )
+		{
+			convertedTaps[ x ] = (float)taps[ x ];
+		}
+		
+		setImpulseTaps( convertedTaps );
+	}
+	
+	public void setImpulseTaps( float[] taps )
+	{
+		double smallest = 0.0;
+		double largest = 0.0;
+		
+        final XYChart.Series<Number, Number> series = new XYChart.Series<>();
+        
+        for( int x = 0; x < taps.length; x++ )
+        {
+        	if( taps[ x ] < smallest )
+			{
+        		smallest = taps[ x ];
+			}
+        	
+        	if( taps[ x ] > largest )
+        	{
+        		largest = taps[ x ];
+        	}
+
+        	series.getData().add( new XYChart.Data<Number, Number>( x, taps[ x ] ) );
+        }
+        
+        mImpulseXAxis.setLowerBound( 0.0 );
+        mImpulseXAxis.setUpperBound( taps.length - 1 );
+        mImpulseXAxis.setTickUnit( 1.0 );
+        
+        double padding = Math.abs( largest - smallest ) * .1;
+        
+        mImpulseYAxis.setLowerBound( smallest - padding );
+        mImpulseYAxis.setUpperBound( largest + padding );
+        mImpulseYAxis.setTickUnit( padding );
+
+        ObservableList<Series<Number,Number>> list = mImpulseChart.getData();
+        list.add( series );
+	}
+	
+	public void setDFTTaps( float[] taps )
+	{
+		double smallest = 0.0;
+		
+        final XYChart.Series<Number, Number> series = new XYChart.Series<>();
+        
+        for( int x = 0; x < taps.length; x++ )
+        {
+        	if( taps[ x ] < smallest )
+			{
+        		smallest = taps[ x ];
+			}
+        	
+        	series.getData().add( new XYChart.Data<Number, Number>( x, taps[ x ] ) );
+        }
+        
+        ObservableList<Series<Number,Number>> list = mDFTChart.getData();
+        list.add( series );
+	}
+	
+    public float[] convertDFTToDecibel( float[] dft, int tapLength )
+    {
+		float[] decibels = new float[ dft.length / 2 ];
+
+		int index = 0;
+		
+		for( int x = 0; x < decibels.length; x ++ )
+		{
+			index = x * 2;
+			
+			decibels[ x ] = 20.0f * (float)Math.log10( ( ( dft[ index ] * dft[ index ] ) + 
+				  ( dft[ index + 1 ] * dft[ index + 1 ] ) ) );
+		}
+
+		return decibels;
+    }
+}

src/dsp/filter/design/FilterViewer.java

@@ -0,0 +1,96 @@
+package dsp.filter.design;
+
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import dsp.filter.fir.FIRFilterSpecification;
+import dsp.filter.fir.remez.RemezFIRFilterDesigner;
+import javafx.application.Application;
+import javafx.scene.Scene;
+import javafx.stage.Stage;
+
+public class FilterViewer extends Application
+{
+	private final static Logger mLog = LoggerFactory.getLogger( FilterViewer.class );
+
+	/**
+	 * Developer tool to visualize filter designs
+	 */
+	public FilterViewer()
+	{
+	}
+	
+	@Override
+	public void start( Stage primaryStage ) throws Exception
+	{
+		Scene scene = new Scene( new FilterView( getFilter() ) );
+		
+		primaryStage.setTitle( "Filter Viewer" );
+		primaryStage.setScene( scene );
+		primaryStage.show();
+	}
+	
+	public static void main( String[] args )
+	{
+		launch( args );
+	}
+
+	/**
+	 * Provides the filter to visualize.  Modify this method to visualize your filter design.
+	 */
+	private float[] getFilter()
+	{
+		FIRFilterSpecification specification = FIRFilterSpecification.lowPassBuilder()
+		        .sampleRate( 10000 )
+		        .gridDensity( 16 )
+                .passBandCutoff( 1300 )
+		        .passBandAmplitude( 1.0 )
+		        .passBandRipple( 0.01 )
+		        .stopBandStart( 2675 )
+		        .stopBandAmplitude( 0.0 )
+		        .stopBandRipple( 0.001 )
+		        .build();
+		
+//		FIRFilterSpecification specification = FIRFilterSpecification.highPassBuilder()
+//		        .sampleRate( 8000 )
+//		        .gridDensity( 16 )
+//		        .stopBandCutoff( 2000 )
+//		        .stopBandAmplitude( 0.0 )
+//		        .stopBandRipple( 0.003 )
+//		        .passBandStart( 3000 )
+//		        .passBandAmplitude( 1.0 )
+//		        .passBandRipple( 0.001 )
+//		        .build();
+
+//		FIRFilterSpecification specification = FIRFilterSpecification.bandPassBuilder()
+//	        .sampleRate( 8000 )
+//	        .order( 164 )
+//	        .stopFrequency1( 1100 )
+//	        .passFrequencyBegin( 1200 )
+//	        .passFrequencyEnd( 2400 )
+//	        .stopFrequency2( 2500 )
+//	        .stopRipple( 0.01 )
+//	        .passRipple( 0.01 )
+//	        .build();
+
+		float[] taps = null;
+
+//		taps = FilterFactory.getLowPass( 16000, 2400, 73, WindowType.HANNING );
+
+		try
+		{
+			RemezFIRFilterDesigner designer = new RemezFIRFilterDesigner( specification );
+			
+			if( designer.isValid() )
+			{
+				taps = designer.getImpulseResponse();
+			}
+		}
+		catch( FilterDesignException fde )
+		{
+		    mLog.error( "Filter design error", fde );
+		}
+		
+		return taps;
+	}
+}

src/dsp/filter/fir/remez/FIRLinearPhaseFilterType.java

@@ -0,0 +1,56 @@
+package dsp.filter.fir.remez;
+
+public enum FIRLinearPhaseFilterType
+{
+	/**
+	 * TYPE 1
+	 *    Length: ODD
+	 *     Order: EVEN
+	 *  Symmetry: SYMMETRICAL
+	 *  
+	 * Valid for: LOW PASS, HIGH PASS
+	 */
+	TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL,
+	
+	/**
+	 * TYPE 2
+	 *    Length: EVEN
+	 *     Order: ODD
+	 *  Symmetry: SYMMETRICAL
+	 * Valid for: LOW PASS
+	 * 
+	 * Notes: 
+	 * -Frequency response at FS/2 (Pi) must be 0 -- not suitable for High Pass.
+	 */
+	TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL,
+	
+	/**
+	 * TYPE 3
+	 *    Length: ODD
+	 *     Order: EVEN
+	 *  Symmetry: ANTI-SYMMETRICAL
+	 * Valid for: 
+	 * 
+	 * Notes:
+	 * -Introduces a 90 degree phase shift such that h[N-n] = -h[n] for all n
+	 * -Frequency response at 0 and FS/2 (Pi) must be 0 -- not suitable for Low or High Pass.
+	 */
+	TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL,
+	
+	/**
+	 * TYPE 4
+	 *    Length: EVEN
+	 *     Order: ODD
+	 *  Symmetry: ANTI-SYMMETRICAL
+	 * Valid for: HIGH PASS
+	 * 
+	 * Notes:
+	 * -Introduces a 90 degree phase shift such that h[N-n] = -h[n] for all n
+	 * -Frequency response at 0 must be zero -- not suitable for Low Pass
+	 */
+	TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL,
+	
+    BANDPASS, 
+    DIFFERENTIATOR, 
+    HILBERT;
+}

src/dsp/filter/fir/remez/Grid.java

@@ -0,0 +1,112 @@
+package dsp.filter.fir.remez;
+
+import dsp.filter.fir.FIRFilterSpecification;
+import dsp.filter.fir.FIRFilterSpecification.FrequencyBand;
+
+public class Grid
+{
+    public static final double TWO_PI = 2.0 * Math.PI;
+    
+    private double[] mFrequencyGrid;
+    private double[] mCosineOfFrequencyGrid;
+    private double[] mDesiredResponse;
+    private double[] mWeight;
+
+    /**
+     * Creates a dense frequency grid from the filter specification for those frequency bands
+     * where the weighted value is non-zero.  This grid is appropriate for use in the Remez
+     * exchange algorithm.
+     */
+    public Grid( FIRFilterSpecification specification )
+    {
+    	create( specification );
+    }
+    
+    private void create( FIRFilterSpecification specification )
+    {
+        int gridSize = specification.getGridSize();
+        
+        mFrequencyGrid = new double[ gridSize ];
+        mCosineOfFrequencyGrid = new double[ gridSize ];
+        mDesiredResponse = new double[ gridSize ];
+        mWeight = new double[ gridSize ];
+        
+        double gridFrequencyInterval = specification.getGridFrequencyInterval();
+        
+        double bandZeroEdge = specification.getFrequencyBands().get( 0 ).getStart();
+        
+        // For differentiator, hilbert, symmetry is odd and Grid[0] = max(delf, bands[0])
+        double grid0 = specification.isAntiSymmetrical() && 
+        			   ( gridFrequencyInterval > bandZeroEdge ) ? gridFrequencyInterval : bandZeroEdge;
+        
+        int j = 0;
+
+        double maxRipple = specification.getMaxBandAmplitude();
+
+        for( int x = 0; x < specification.getFrequencyBands().size(); x++ )
+        {
+        	FrequencyBand band = specification.getFrequencyBands().get( x );
+        	
+            double lowFrequency = ( x == 0 ? grid0 : band.getStart() );
+
+            for( int i = 0; i < band.getGridSize(); i++ )
+        	{
+                mDesiredResponse[ j ] = band.getAmplitude();
+
+                //Set the weighting function equally across this band's dense grids
+                mWeight[ j ] = band.getWeight( maxRipple );
+
+                //Frequency grid is linear taper from start frequency to stop frequency in
+                //gridFrequencyInterval increments
+                mFrequencyGrid[ j ] = lowFrequency;
+                
+                mCosineOfFrequencyGrid[ j ] = Math.cos( TWO_PI * mFrequencyGrid[ j ] );
+
+                lowFrequency += gridFrequencyInterval;
+                
+                j++;
+        	}
+            
+            int lastIndex = j - 1;
+            
+            mFrequencyGrid[ lastIndex ] = band.getEnd();
+        }
+    }
+
+    /**
+     * Size of the dense frequency grid
+     */
+    public int getSize()
+    {
+        return mFrequencyGrid.length;
+    }
+
+    /**
+     * The dense frequency grid
+     */
+    public double[] getFrequencyGrid()
+    {
+        return mFrequencyGrid;
+    }
+    
+    public double[] getCosineFrequencyGrid()
+    {
+    	return mCosineOfFrequencyGrid;
+    }
+
+    /**
+     * Desired frequency response
+     */
+    public double[] getDesiredResponse()
+    {
+        return mDesiredResponse;
+    }
+
+    /**
+     * Weighting function
+     */
+    public double[] getWeight()
+    {
+        return mWeight;
+    }
+}

src/dsp/filter/fir/remez/RemezFIRFilterDesigner.java

@@ -0,0 +1,675 @@
+package dsp.filter.fir.remez;
+
+import java.text.DecimalFormat;
+import java.util.ArrayList;
+import java.util.Iterator;
+import java.util.List;
+
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import dsp.filter.design.FilterDesignException;
+import dsp.filter.fir.FIRFilterSpecification;
+
+public class RemezFIRFilterDesigner
+{
+	private final static Logger mLog = LoggerFactory.getLogger( RemezFIRFilterDesigner.class );
+	
+	private static final double CONVERGENCE_THRESHOLD = 0.0001;
+    public static final int MAXIMUM_ITERATION_COUNT = 40;
+    public static final double TWO_PI = 2.0 * Math.PI;
+
+    private static final DecimalFormat DECIMAL_FORMATTER = new DecimalFormat( "0.00000000" );
+
+    private FIRFilterSpecification mSpecification;
+	private Grid mGrid;
+	private List<Integer> mExtremalIndices;
+	
+	private double[] mD;
+	private double[] mGridErrors;
+	private double[] mGridFrequencyResponse;
+	private double[] mIdealFrequencyResponse;
+	
+	private double mDelta;
+	
+	private boolean mConverged;
+
+	/**
+	 * Linear FIR filter designer using the Parks-McClelland Remez exchange algorithm.
+	 * 
+	 * Construct a filter specification using the builder methods in the FIRFilterSpecification 
+	 * class and construct an instance of this class with the specification.  Check the isValid()
+	 * method after construction to see if the filter design was successful.  If successful, use the
+	 * getFrequencyResponse() and getImpulseResponse() methods to access the filter taps as the 
+	 * specified length.
+	 * 
+	 * If the filter cannot be designed from the specification the isValid() method will return 
+	 * false and any attempts to access the getFrequencyResponse() or getImpulseResponse() methods
+	 * will throw a FilterDesignException.
+	 * 
+	 * @param specification that defines a Type 1-4 linear phase FIR Filter
+	 */
+	public RemezFIRFilterDesigner( FIRFilterSpecification specification )
+	{
+		mSpecification = specification;
+		
+		design();
+	}
+
+	/**
+	 * Designs the filter according to the specification
+	 */
+	private void design()
+	{
+		mGrid = new Grid( mSpecification );
+		mExtremalIndices = getInitialExtremalIndices();
+
+    	int iterationCount = 0;
+    	
+    	try
+    	{
+        	do
+        	{
+            	calculateGridFrequencyResponse();
+            	calculateGridError();
+            	findExtremalIndices();
+            	checkConvergence();
+            	
+            	iterationCount++;
+        	}
+        	while( !mConverged && iterationCount < MAXIMUM_ITERATION_COUNT );
+    	}
+    	catch( FilterDesignException fde )
+    	{
+    		mLog.error( "Filter design error - couldn't find extremal indices at count [" + 
+    				iterationCount + "] try changing filter order up/down from [" + 
+    				mSpecification.getOrder() + "]" );
+    		
+    		mConverged = false;
+    	}
+
+    	if( mConverged )
+    	{
+        	//Update frequency response using the final set of extremal indices
+        	calculateGridFrequencyResponse();
+    	}
+	}
+
+	/**
+	 * Indicates if the filter was successfully created from the filter specification.  Check this
+	 * method before accessing frequency response or filter response, otherwise a FilterDesignException
+	 * is thrown.
+	 * 
+	 * @return true if the filter was successfully designed.
+	 */
+	public boolean isValid()
+	{
+		return mConverged;
+	}
+
+	/**
+	 * Frequency response of the designed filter with double precision
+	 * 
+	 * @return half ( length / 2 ) of filter frequency response sampled at 0.0 - PI radians
+	 * @throws FilterDesignException if the specified filter cannot be designed
+	 */
+	public double[] getFrequencyResponse() throws FilterDesignException
+	{
+		if( !mConverged )
+		{
+			throw new FilterDesignException( "Can't create filter from specification - failed "
+					+ "to converge" );
+		}
+		
+		//Resample the final polynomial at the desired filter length
+    	double[] resampled = resample();
+
+    	//Apply Parks/McClelland frequency response corrections according to the filter type
+    	return correctFrequencyResponse( resampled );
+	}
+
+	/**
+	 * Impulse response of the designed filter with floating point precision.
+	 * 
+	 * @return filter impulse response
+	 * @throws FilterDesignException if the specified filter cannot be designed
+	 */
+    public float[] getImpulseResponse() throws FilterDesignException
+	{
+		return convertToFloatArray( getImpulseResponseDoubles() );
+	}
+	
+    /**
+     * Impulse response of the designed filter with double precision
+     */
+    public double[] getImpulseResponseDoubles() throws FilterDesignException
+    {
+    	double[] frequencyResponse = getFrequencyResponse();
+    	
+    	int length = mSpecification.getFilterLength();
+    	
+    	double[] impulseResponse = new double[ length ];
+
+    	double M;
+    	
+    	switch( mSpecification.getFilterType() )
+    	{
+    		case TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL:
+    			M = ( (double)length - 1.0 ) / 2.0;
+    	    	
+    			for( int n = 0; n < length; n++ )
+    			{
+    				double accumulator = frequencyResponse[ 0 ];
+    				
+    				double frequency = TWO_PI * ( n - M ) / length;
+    				
+    				for( int k = 1; k <= M; k++ )
+    				{
+    					accumulator += 2.0 * frequencyResponse[ k ] * Math.cos( frequency * (double)k );
+    				}
+    				
+    				impulseResponse[ n ] = accumulator / (double)length;
+    			}
+    			break;
+    		case TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL:
+    			double offset = (double)( length - 1 ) / 2.0;
+    			
+                for( int n = 0; n < length; n++ )
+                {
+    				double accumulator = frequencyResponse[ 0 ];
+                    
+    				double frequency = TWO_PI * ( (double)n - offset ) / (double)length;
+
+                    for( int k = 1; k < frequencyResponse.length; k++ )
+                    {
+                    	accumulator += 2.0 * frequencyResponse[ k ] * Math.cos( frequency * (double)k );
+                    }
+
+                    impulseResponse[ n ] = accumulator / (double)length;
+                }
+                break;
+    	}
+    	
+    	return impulseResponse;
+    }
+
+    /**
+     * Calculates the frequency response of this filter at the specified cosine of the frequency
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 116a to
+     * determine the actual response of the specified frequency cosine using Lagrange 
+     * interpolation of the polynomial formed at the extremal index set.
+     * 
+     * @param cosineOfFrequency to calculate the frequency response
+     * 
+     * @return actual frequency response for the specified frequency with the current set of
+     * extreaml indices
+     */
+    private double getFrequencyResponse( double cosineOfFrequency )
+    {
+    	double numerator = 0.0;
+    	double denominator = 0.0;
+    	
+		for( int k = 0; k < mExtremalIndices.size() - 1; k++ )
+    	{
+    		double cosineDelta = cosineOfFrequency - 
+    				mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( k ) ];
+    		
+    		//If this frequency is or is close to one of the polynomial points, use c[k] for that point
+    		if( Math.abs( cosineDelta ) < 1.0e-7 )
+    		{
+    			numerator = mIdealFrequencyResponse[ k ];
+    			denominator = 1.0;
+    			break;
+    		}
+    		else
+    		{
+    			double dkOverCosineDelta = mD[ k ] / cosineDelta;
+    			
+        		numerator += dkOverCosineDelta * mIdealFrequencyResponse[ k ];
+        		denominator += dkOverCosineDelta;
+    		}
+    	}
+
+		return numerator / denominator;
+    }
+    
+    /**
+     * Calculated Delta value 
+     * 
+     * @return delta value
+     */
+	public double getDelta()
+	{
+		return mDelta;
+	}
+	
+    /**
+     * Creates initial array of extremal frequency indices spaced at intervals of the grid frequency
+     * set for all frequencies where a frequency band weighting is specified.
+     */
+    private List<Integer> getInitialExtremalIndices()
+    {
+    	int count = mSpecification.getExtremaCount();
+    	
+    	int density = mSpecification.getGridDensity();
+    	
+    	List<Integer> extremalIndices = new ArrayList<>();
+
+        for( int i = 0; i < count; i++ )
+        {
+            extremalIndices.add( i * density );
+        }
+
+        return extremalIndices;
+    }
+    
+    /**
+     * Calculates the actual frequency response of a dense frequency grid using the polynomial 
+     * represented by the extremal indices.
+     */
+    private void calculateGridFrequencyResponse()
+    {
+    	double[] b = calculateB();
+    	
+    	calculateDelta( b );
+    	
+    	calculateC();
+    	
+    	calculateD( b );
+
+    	updateGridFrequencyResponse();
+    }
+
+    /**
+     * Updates the frequency response of frequencies specified in the dense frequency grid
+     */
+    private void updateGridFrequencyResponse()
+    {
+    	mGridFrequencyResponse = new double[ mGrid.getSize() ];
+
+    	double[] gridFrequencyCosines = mGrid.getCosineFrequencyGrid();
+    	
+    	for( int i = 0; i < mGridFrequencyResponse.length; i++ )
+    	{
+    		mGridFrequencyResponse[ i ] = getFrequencyResponse( gridFrequencyCosines[ i ] );
+    	}
+    }
+	
+    /**
+     * Calculates the value of b across the set of extremal indices.
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 115
+     * 
+     * @return the array of b values of length L+2
+     */
+    private double[] calculateB()
+    {
+    	int length = mExtremalIndices.size();
+    	
+    	double[] b = new double[ length ];
+
+    	for( int k = 0; k < length; k++ )
+    	{
+    		b[ k ] = 1.0;
+    		
+    		double xk = mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( k ) ];
+    		
+        	for( int i = 0; i < length; i++ )
+        	{
+        		if( i != k )
+        		{
+        			double xi = mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( i ) ];
+        			
+        			double denominator = xk - xi;
+        			
+                    if( Math.abs( denominator ) < 0.00001 )
+                    {
+                        denominator = 0.00001;
+                    }
+                    
+        			b[ k ] *= 1.0 / denominator;
+        		}
+        	}
+    	}
+    	
+    	return b;
+    }
+    
+    /**
+     * Calculates the value of delta that represents the maximum ripple for the current set of 
+     * extremal indices (L+2)
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 114
+     * 
+     * @param b calculated for the grid from equation 115
+     */
+    private void calculateDelta( double[] b )
+    {
+    	double numerator = 0.0;
+    	double denominator = 0.0;
+
+    	double sign = 1.0;
+    	
+    	for( int k = 0; k < b.length; k++ )
+    	{
+    		if( k < mExtremalIndices.size() )
+    		{
+    			int extremalIndex = mExtremalIndices.get( k );
+    			
+        		numerator += ( b[ k ] * mGrid.getDesiredResponse()[ extremalIndex ] );
+        		denominator += b[ k ] * sign / mGrid.getWeight()[ extremalIndex ];
+        		sign = -sign;
+    		}
+    		else
+    		{
+    			mLog.error( "Something went wrong -- the length of b exceeds the set of extremal indices" );
+    		}
+    	}
+    	
+    	mDelta = numerator / denominator;
+    }
+
+    /**
+     * Calculates the ideal frequency response (C) of the current extremal index set.  This reference
+     * is used to determine the level of error in the set of actual frequency responses calculated
+     * from this set of extremal indices to determine convergence.
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 116b
+     */
+    private void calculateC()
+    {
+    	int length = mSpecification.getExtremaCount() - 1;
+    	
+    	mIdealFrequencyResponse = new double[ length ];
+
+    	double sign = 1.0;
+    	
+    	for( int k = 0; k < length; k++ )
+    	{
+    		if( k < mExtremalIndices.size() )
+    		{
+    			int index = mExtremalIndices.get( k );
+    			
+        		mIdealFrequencyResponse[ k ] = mGrid.getDesiredResponse()[ index ] - 
+        				( sign * mDelta / mGrid.getWeight()[ index ] );
+        		
+        		sign = -sign;
+    		}
+    	}
+    }
+
+    /**
+     * Calculates the set of D values for the current extremal index set.
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 116c
+     */
+    private void calculateD( double[] b )
+    {
+    	int length = mExtremalIndices.size() - 1;
+    	
+    	mD = new double[ length ];
+    	
+    	for( int k = 0; k < length; k++ )
+    	{
+        	//Note: extremalCosines array is one index longer than d array, so we use length as
+    		//a pointer to the final (L+2) extremalCosine index
+    		mD[k] = b[k] * ( mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( k ) ] - 
+    						 mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( length ) ] );
+    	}
+    }
+
+    /**
+     * Calculates the weighted error between the specification desired frequency response and the 
+     * actual frequency response of the current extremal index set across the frequencies in the
+     * dense frequency grid.
+     * 
+     * Implements Oppenheim/Schafer Discrete Time Signal Processing, 3e, 2016, equation 112
+     */
+    public void calculateGridError()
+    {
+    	int length = mGridFrequencyResponse.length;
+    	
+    	mGridErrors = new double[ length ];
+    	
+    	for( int i = 0; i < length; i++ )
+    	{
+            mGridErrors[ i ] = mGrid.getWeight()[ i ] * 
+            		( mGrid.getDesiredResponse()[ i ] - mGridFrequencyResponse[ i ] );
+    	}
+    }
+
+    /**
+     * Identifies the indexes of the extremal error values using the Alternation theorem
+     * 
+     * @throws FilterDesignException if the full set of extremal indices can't be identified
+     */
+    private void findExtremalIndices() throws FilterDesignException
+    {
+    	mExtremalIndices.clear();
+    	
+    	//Check for extremum at error index 0
+    	if( ( ( mGridErrors[ 0 ] > 0.0 && mGridErrors[ 0 ] > mGridErrors[ 1 ] ) ||
+    		  ( mGridErrors[ 0 ] < 0.0 && mGridErrors[ 0 ] < mGridErrors[ 1 ] ) ) &&
+    			isGTEDelta( mGridErrors[ 0 ] ) )
+		{
+    		mExtremalIndices.add( 0 );
+		}
+    	
+        //Check for extrema in middle error indices )
+        for( int x = 1; x < mGridErrors.length - 1; x++ )
+        {
+        	if( ( ( mGridErrors[ x ] > 0.0 && ( mGridErrors[ x - 1 ] <= mGridErrors[ x ] && mGridErrors[ x ] > mGridErrors[ x + 1 ] ) ) ||
+        		  ( mGridErrors[ x ] < 0.0 && ( mGridErrors[ x - 1 ] >= mGridErrors[ x ] && mGridErrors[ x ] < mGridErrors[ x + 1 ] ) ) ) &&
+        			isGTEDelta( mGridErrors[ x ] ) )
+    		{
+        		mExtremalIndices.add( x );
+    		}
+        }
+        
+        //Check for extremum at final error index
+        int last = mGridErrors.length - 1;
+        
+        if( ( ( mGridErrors[ last ] > 0.0 && ( mGridErrors[ last ] > mGridErrors[ last - 1 ] ) ) ||
+        	  ( mGridErrors[ last ] < 0.0 && ( mGridErrors[ last ] < mGridErrors[ last - 1 ] ) ) ) &&
+        		isGTEDelta( mGridErrors[ last ] ) )
+        {
+        	mExtremalIndices.add( last );
+        }
+
+        //Ensure we have the minimum of extremals before we continue
+        if( mExtremalIndices.size() < mSpecification.getExtremaCount() )
+        {
+            throw new FilterDesignException( "Couldn't find the minimum extremal frequencies in "
+            		+ "error set" );
+        }
+
+        //Enforce alternation theory -- only one extremal for each transition about the zero axis 
+        List<Integer> indicesToRemove = new ArrayList<>();
+        
+        Iterator<Integer> it = mExtremalIndices.iterator();
+        Integer current = it.next();
+        Integer next;
+
+        boolean positiveAxis = mGridErrors[ current ] > 0.0;
+
+        while( it.hasNext() )
+        {
+        	next = it.next();
+        	
+        	//Check for consecutive (redundant) indices on same side of zero axis - retain largest
+        	if( !( positiveAxis ^ ( mGridErrors[ next ] > 0.0 ) ) )
+        	{
+        		if( Math.abs( mGridErrors[ next ] ) <= Math.abs( mGridErrors[ current ] ))
+        		{
+        			//Remove next error index that is less than current index and on 
+        			it.remove();
+        			next = current;
+        		}
+        		else
+        		{
+        			//Since we can't remove the current index with the iterator, queue for later removal
+        			indicesToRemove.add( current );
+        		}
+        	}
+        	else
+        	{
+        		positiveAxis = !positiveAxis;
+        	}
+        	
+        	current = next;
+        }
+
+        //Remove redundant extremal indices that couldn't be removed via the iterator
+    	mExtremalIndices.removeAll( indicesToRemove );
+
+        //Truncate the list to one larger than needed by removing excess tailing indices
+        while( mExtremalIndices.size() > mSpecification.getExtremaCount() )
+        {
+        	mExtremalIndices.remove( mExtremalIndices.size() - 1 );
+        }
+        
+        //If we have one too many indices, delete the smaller of the first or last
+        if( mExtremalIndices.size() > mSpecification.getExtremaCount() )
+        {
+        	int lastIndex = mExtremalIndices.size() - 1;
+        	
+        	if( Math.abs( mGridErrors[ mExtremalIndices.get( 0 )] ) > 
+        		Math.abs( mGridErrors[ mExtremalIndices.get( lastIndex )] ) )
+			{
+        		mExtremalIndices.remove( lastIndex );
+			}
+        	else
+        	{
+        		mExtremalIndices.remove( 0 );
+        	}
+        }
+        
+        //Detect if we have too few indices
+        if( mExtremalIndices.size() < mSpecification.getExtremaCount() )
+        {
+            throw new FilterDesignException( "Couldn't find the minimum extremal frequencies in "
+            		+ "error set" );
+        }
+    }
+
+    /**
+     * Indicates if the absolute value of the argument is greater than or equal to the delta
+     * value with accuracy to 14 digits of precision.  This avoids rounding errors at 15 digits or 
+     * precision or greater.
+     * 
+     * @param value to evaluate
+     * 
+     * @return true if the absolute value is greater than or equal to the absolute delta value
+     */
+    private boolean isGTEDelta( double value )
+    {
+    	return Math.abs( value ) - Math.abs( mDelta ) > -1.0e-5;
+    }
+    
+    /**
+     * Checks for convergence of the frequency response of the current set of extremal indices to
+     * the filter specification by comparing the maximum absolute error value against the delta
+     * value to determine if these two values are within the convergence threshold.
+     */
+    public void checkConvergence()
+    {
+        double maximum = Math.abs( mGridErrors[ mExtremalIndices.get( 0 ) ] );
+
+        for( int i = 1; i < mExtremalIndices.size(); i++ )
+        {
+            double current = Math.abs( mGridErrors[ mExtremalIndices.get( i ) ] );
+
+            if( current > maximum )
+            {
+                maximum = current;
+            }
+        }
+        
+        double convergence = maximum - Math.abs( mDelta );
+        
+        mConverged = convergence < CONVERGENCE_THRESHOLD;
+    }
+    
+    /**
+     * Coverts/casts the double array to a float array 
+     */
+    private static float[] convertToFloatArray( double[] samples )
+    {
+    	float[] converted = new float[ samples.length ];
+    	
+    	for( int x = 0; x < samples.length; x++ )
+    	{
+    		converted[ x ] = (float)samples[ x ];
+    	}
+    	
+    	return converted;
+    }
+    
+    /**
+     * Resamples the designed filter at half of the filter length at evenly spaced intervals of the
+     * frequency spectrum from 0 to PI radians. 
+     * 
+     * @return frequency response of the filter sampled the desired filter length and order
+     */
+    private double[] resample()
+    {
+    	int length = mSpecification.getFilterLength();
+
+    	if( length % 2 == 0 )
+    	{
+    		length--;
+    	}
+
+    	double half = (double)length / 2.0;
+
+    	double[] resampled = new double[ (int)Math.ceil( half ) ] ;
+    	
+    	for( int x = 0; x < resampled.length; x++ )
+    	{
+    		resampled[ x ] = getFrequencyResponse( Math.cos( Math.PI * (double)x / half ) );
+    	}
+
+    	return resampled;
+    }
+    
+    /**
+     * Applies Parks/McClelland filter type correction to the final frequency response array
+     * 
+     * @param frequencyResponse of the resampled filter design from 0 <> ~ PI radians
+     * 
+     * @return corrected frequency response
+     */
+    private double[] correctFrequencyResponse( double[] frequencyResponse )
+    {
+    	double filterLength = mSpecification.getFilterLength();
+    	
+    	switch( mSpecification.getFilterType() )
+    	{
+			case TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL:
+				//No correction needed
+				break;
+			case TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL:
+				//No correction needed
+				break;
+			case TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL:
+				for( int x = 0; x < frequencyResponse.length; x++ )
+				{
+		    		double frequencyRadians = Math.PI * ( ( (double)x ) / filterLength );
+					frequencyResponse[ x ] *= Math.sin( TWO_PI * frequencyRadians );
+				}
+				break;
+			case TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL:
+				for( int x = 0; x < frequencyResponse.length; x++ )
+				{
+		    		double frequencyRadians = Math.PI * ( ( (double)x ) / filterLength );
+					frequencyResponse[ x ] *= Math.sin( Math.PI * frequencyRadians );
+				}
+				break;
+			default:
+				break;
+    	}
+    	
+    	return frequencyResponse;
+    }
+}

src/gui/control/CurveFittedAreaChart.java

@@ -0,0 +1,207 @@
+package gui.control;
+
+import javafx.collections.ObservableList;
+import javafx.geometry.Point2D;
+import javafx.scene.Group;
+import javafx.scene.chart.AreaChart;
+import javafx.scene.chart.NumberAxis;
+import javafx.scene.chart.XYChart;
+import javafx.scene.shape.ClosePath;
+import javafx.scene.shape.CubicCurveTo;
+import javafx.scene.shape.LineTo;
+import javafx.scene.shape.MoveTo;
+import javafx.scene.shape.Path;
+import javafx.scene.shape.PathElement;
+import javafx.util.Pair;
+
+public class CurveFittedAreaChart extends AreaChart<Number, Number>
+{
+
+	public CurveFittedAreaChart( NumberAxis xAxis, NumberAxis yAxis )
+	{
+		super( xAxis, yAxis );
+	}
+
+	@Override
+	protected void layoutPlotChildren()
+	{
+		super.layoutPlotChildren();
+		for ( int seriesIndex = 0; seriesIndex < getDataSize(); seriesIndex++ )
+		{
+			final XYChart.Series<Number, Number> series = getData()
+					.get( seriesIndex );
+			final Path seriesLine = (Path) ( (Group) series.getNode() )
+					.getChildren().get( 1 );
+			final Path fillPath = (Path) ( (Group) series.getNode() )
+					.getChildren().get( 0 );
+			smooth( seriesLine.getElements(), fillPath.getElements() );
+		}
+	}
+
+	private int getDataSize()
+	{
+		final ObservableList<XYChart.Series<Number, Number>> data = getData();
+		return ( data != null ) ? data.size() : 0;
+	}
+
+	private static void smooth( ObservableList<PathElement> strokeElements,
+			ObservableList<PathElement> fillElements )
+	{
+		// as we do not have direct access to the data, first recreate the list
+		// of all the data points we have
+		final Point2D[] dataPoints = new Point2D[ strokeElements.size() ];
+		for ( int i = 0; i < strokeElements.size(); i++ )
+		{
+			final PathElement element = strokeElements.get( i );
+			if ( element instanceof MoveTo )
+			{
+				final MoveTo move = (MoveTo) element;
+				dataPoints[ i ] = new Point2D( move.getX(), move.getY() );
+			} else if ( element instanceof LineTo )
+			{
+				final LineTo line = (LineTo) element;
+				final double x = line.getX(), y = line.getY();
+				dataPoints[ i ] = new Point2D( x, y );
+			}
+		}
+		// next we need to know the zero Y value
+		final double zeroY = ( (MoveTo) fillElements.get( 0 ) ).getY();
+		// now clear and rebuild elements
+		strokeElements.clear();
+		fillElements.clear();
+		Pair<Point2D[], Point2D[]> result = calcCurveControlPoints(
+				dataPoints );
+		Point2D[] firstControlPoints = result.getKey();
+		Point2D[] secondControlPoints = result.getValue();
+		// start both paths
+		strokeElements.add(
+				new MoveTo( dataPoints[ 0 ].getX(), dataPoints[ 0 ].getY() ) );
+		fillElements.add( new MoveTo( dataPoints[ 0 ].getX(), zeroY ) );
+		fillElements.add(
+				new LineTo( dataPoints[ 0 ].getX(), dataPoints[ 0 ].getY() ) );
+		// add curves
+		for ( int i = 1; i < dataPoints.length; i++ )
+		{
+			final int ci = i - 1;
+			strokeElements
+					.add( new CubicCurveTo( firstControlPoints[ ci ].getX(),
+							firstControlPoints[ ci ].getY(),
+							secondControlPoints[ ci ].getX(),
+							secondControlPoints[ ci ].getY(),
+							dataPoints[ i ].getX(), dataPoints[ i ].getY() ) );
+			fillElements.add( new CubicCurveTo( firstControlPoints[ ci ].getX(),
+					firstControlPoints[ ci ].getY(),
+					secondControlPoints[ ci ].getX(),
+					secondControlPoints[ ci ].getY(), dataPoints[ i ].getX(),
+					dataPoints[ i ].getY() ) );
+		}
+		// end the paths
+		fillElements.add( new LineTo(
+				dataPoints[ dataPoints.length - 1 ].getX(), zeroY ) );
+		fillElements.add( new ClosePath() );
+	}
+
+	/**
+	 * Calculate open-ended Bezier Spline Control Points.
+	 */
+	public static Pair<Point2D[], Point2D[]> calcCurveControlPoints(
+			Point2D[] dataPoints )
+	{
+		Point2D[] firstControlPoints;
+		Point2D[] secondControlPoints;
+		int n = dataPoints.length - 1;
+		if ( n == 1 )
+		{ // Special case: Bezier curve should be a straight line.
+			firstControlPoints = new Point2D[ 1 ];
+			// 3P1 = 2P0 + P3
+			firstControlPoints[ 0 ] = new Point2D(
+					( 2 * dataPoints[ 0 ].getX() + dataPoints[ 1 ].getX() ) / 3,
+					( 2 * dataPoints[ 0 ].getY() + dataPoints[ 1 ].getY() )
+							/ 3 );
+
+			secondControlPoints = new Point2D[ 1 ];
+			// P2 = 2P1 – P0
+			secondControlPoints[ 0 ] = new Point2D(
+					2 * firstControlPoints[ 0 ].getX() - dataPoints[ 0 ].getX(),
+					2 * firstControlPoints[ 0 ].getY()
+							- dataPoints[ 0 ].getY() );
+			return new Pair<Point2D[], Point2D[]>( firstControlPoints,
+					secondControlPoints );
+		}
+
+		// Calculate first Bezier control points
+		// Right hand side vector
+		double[] rhs = new double[ n ];
+
+		// Set right hand side X values
+		for ( int i = 1; i < n - 1; ++i )
+		{
+			rhs[ i ] = 4 * dataPoints[ i ].getX()
+					+ 2 * dataPoints[ i + 1 ].getX();
+		}
+		rhs[ 0 ] = dataPoints[ 0 ].getX() + 2 * dataPoints[ 1 ].getX();
+		rhs[ n - 1 ] = ( 8 * dataPoints[ n - 1 ].getX()
+				+ dataPoints[ n ].getX() ) / 2.0;
+		// Get first control points X-values
+		double[] x = GetFirstControlPoints( rhs );
+
+		// Set right hand side Y values
+		for ( int i = 1; i < n - 1; ++i )
+		{
+			rhs[ i ] = 4 * dataPoints[ i ].getY()
+					+ 2 * dataPoints[ i + 1 ].getY();
+		}
+		rhs[ 0 ] = dataPoints[ 0 ].getY() + 2 * dataPoints[ 1 ].getY();
+		rhs[ n - 1 ] = ( 8 * dataPoints[ n - 1 ].getY()
+				+ dataPoints[ n ].getY() ) / 2.0;
+		// Get first control points Y-values
+		double[] y = GetFirstControlPoints( rhs );
+
+		// Fill output arrays.
+		firstControlPoints = new Point2D[ n ];
+		secondControlPoints = new Point2D[ n ];
+		for ( int i = 0; i < n; ++i )
+		{
+			// First control point
+			firstControlPoints[ i ] = new Point2D( x[ i ], y[ i ] );
+			// Second control point
+			if ( i < n - 1 )
+			{
+				secondControlPoints[ i ] = new Point2D(
+						2 * dataPoints[ i + 1 ].getX() - x[ i + 1 ],
+						2 * dataPoints[ i + 1 ].getY() - y[ i + 1 ] );
+			} else
+			{
+				secondControlPoints[ i ] = new Point2D(
+						( dataPoints[ n ].getX() + x[ n - 1 ] ) / 2,
+						( dataPoints[ n ].getY() + y[ n - 1 ] ) / 2 );
+			}
+		}
+		return new Pair<Point2D[], Point2D[]>( firstControlPoints,
+				secondControlPoints );
+	}
+
+	/**
+	 * Solves a tridiagonal system for one of coordinates (x or y) of first
+	 * Bezier control points.
+	 */
+	private static double[] GetFirstControlPoints( double[] rhs )
+	{
+		int n = rhs.length;
+		double[] x = new double[ n ]; // Solution vector.
+		double[] tmp = new double[ n ]; // Temp workspace.
+		double b = 2.0;
+		x[ 0 ] = rhs[ 0 ] / b;
+		for ( int i = 1; i < n; i++ )
+		{// Decomposition and forward substitution.
+			tmp[ i ] = 1 / b;
+			b = ( i < n - 1 ? 4.0 : 3.5 ) - tmp[ i ];
+			x[ i ] = ( rhs[ i ] - x[ i - 1 ] ) / b;
+		}
+		for ( int i = 1; i < n; i++ )
+		{
+			x[ n - i - 1 ] -= tmp[ n - i ] * x[ n - i ]; // Backsubstitution.
+		}
+		return x;
+	}
+}