diff --git a/src/dsp/filter/design/FilterViewer.java b/src/dsp/filter/design/FilterViewer.java
index f5d70995..348645b3 100644
--- a/src/dsp/filter/design/FilterViewer.java
+++ b/src/dsp/filter/design/FilterViewer.java
@@ -40,26 +40,26 @@ public class FilterViewer extends Application
*/
private float[] getFilter()
{
- FIRFilterSpecification specification = FIRFilterSpecification.lowPassBuilder()
- .sampleRate(100000)
- .gridDensity(16)
- .passBandAmplitude(1.0)
- .passBandCutoff(12000)
- .passBandRipple(0.01)
- .stopBandAmplitude(0.0)
- .stopBandRipple(0.001)
- .stopBandStart(13000)
- .build();
+// FIRFilterSpecification specification = FIRFilterSpecification.lowPassBuilder()
+// .sampleRate(10000)
+// .gridDensity(16)
+// .passBandAmplitude(1.0)
+// .passBandCutoff(4000)
+// .passBandRipple(0.01)
+// .stopBandAmplitude(0.0)
+// .stopBandRipple(0.001)
+// .stopBandStart(5000)
+// .build();
-// FIRFilterSpecification specification = FIRFilterSpecification.highPassBuilder()
-// .sampleRate( 24000 )
-// .stopBandCutoff( 800 )
-// .stopBandAmplitude( 0.0 )
-// .stopBandRipple( 0.03 )
-// .passBandStart( 1000 )
-// .passBandAmplitude( 1.0 )
-// .passBandRipple( 0.08 )
-// .build();
+ FIRFilterSpecification specification = FIRFilterSpecification.highPassBuilder()
+ .sampleRate( 24000 )
+ .stopBandCutoff( 800 )
+ .stopBandAmplitude( 0.0 )
+ .stopBandRipple( 0.03 )
+ .passBandStart( 1000 )
+ .passBandAmplitude( 1.0 )
+ .passBandRipple( 0.08 )
+ .build();
// FIRFilterSpecification specification = FIRFilterSpecification.bandPassBuilder()
// .sampleRate( 48000 )
diff --git a/src/dsp/filter/fir/FIRFilterSpecification.java b/src/dsp/filter/fir/FIRFilterSpecification.java
index b03e82ec..ab599f89 100644
--- a/src/dsp/filter/fir/FIRFilterSpecification.java
+++ b/src/dsp/filter/fir/FIRFilterSpecification.java
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * sdrtrunk
+ * Copyright (C) 2014-2017 Dennis Sheirer
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see
+ *
+ ******************************************************************************/
package dsp.filter.fir;
import java.text.DecimalFormat;
@@ -19,14 +37,12 @@ public class FIRFilterSpecification
/**
* Constructs a Remez filter specification. Use one of the filter builders to construct a Remez
* filter specification.
- *
+ *
+ * @param type of filter
* @param order of the filter
- * @param frequencyBands containing the start and stop values for each frequency band normalized
- * to the sample rate where each frequency value is in the range (0.0 <> 0.5)
- * @param amplitudes of each frequency band
- * @param weights of the band ripples normalized to the largest ripple value
+ * @param gridDensity to use when calculating the coefficients
*/
- private FIRFilterSpecification( FIRLinearPhaseFilterType type, int order, int gridDensity )
+ protected FIRFilterSpecification( FIRLinearPhaseFilterType type, int order, int gridDensity )
{
mRemezFilterType = type;
mOrder = order;
@@ -70,6 +86,14 @@ public class FIRFilterSpecification
updateGridSize();
}
+ /**
+ * Removes all frequency bands from this filter specification
+ */
+ public void clearFrequencyBands()
+ {
+ mFrequencyBands.clear();
+ }
+
public List getFrequencyBands()
{
return mFrequencyBands;
@@ -868,7 +892,139 @@ public class FIRFilterSpecification
}
}
- // /**
+ /**
+ * Builder for creating polyphase channelizer low-pass prototype filter.
+ */
+ public static class PolyphaseChannelizerBuilder
+ {
+ private int mOrder;
+ private int mGridDensity = 16;
+ private int mSampleRate;
+ private int mChannelBandwidth;
+ private int mPassFrequencyBegin = 0;
+ private int mPassFrequencyEnd;
+ private int mStopFrequencyBegin;
+ private double mStopRipple = 0.001;
+ private double mPassRipple = 0.01;
+ private double mStopAmplitude = 0.0;
+ private double mPassAmplitude = 1.0;
+
+ public PolyphaseChannelizerBuilder()
+ {
+ }
+
+ /**
+ * Sets the filter sample rate in hertz
+ */
+ public PolyphaseChannelizerBuilder sampleRate(int sampleRate )
+ {
+ mSampleRate = sampleRate;
+ return this;
+ }
+
+ /**
+ * Sets the channel bandwidth in hertz
+ */
+ public PolyphaseChannelizerBuilder channelBandwidth(int channelBandwidth)
+ {
+ mChannelBandwidth = channelBandwidth;
+ return this;
+ }
+
+ /**
+ * Sets the filter order. If the order is not specified, it will be calculated from the
+ * other parameters. If you specify and odd-length order, it will be corrected (+1) to make
+ * it an even order number.
+ */
+ public PolyphaseChannelizerBuilder order(int order )
+ {
+ mOrder = order;
+ return this;
+ }
+
+ /**
+ * Sets the frequency grid density. Default is 16.
+ */
+ public PolyphaseChannelizerBuilder gridDensity(int density )
+ {
+ mGridDensity = density;
+ return this;
+ }
+
+ /**
+ * Specifies the ripple in dB for both stop bands
+ */
+ public PolyphaseChannelizerBuilder stopRipple(double rippleDb )
+ {
+ mStopRipple = rippleDb;
+ return this;
+ }
+
+ /**
+ * Specifies the ripple in dB for the pass band
+ */
+ public PolyphaseChannelizerBuilder passRipple(double rippleDb )
+ {
+ mPassRipple = rippleDb;
+ return this;
+ }
+
+ /**
+ * Specifies the amplitude of the stop bands. Default is 0.0
+ */
+ public PolyphaseChannelizerBuilder stopAmplitude(double amplitude )
+ {
+ mStopAmplitude = amplitude;
+ return this;
+ }
+
+ /**
+ * Specifies the amplitude of the pass band. Default is 1.0
+ */
+ public PolyphaseChannelizerBuilder passAmplitude(double amplitude )
+ {
+ mPassAmplitude = amplitude;
+ return this;
+ }
+
+ public FIRFilterSpecification build()
+ {
+ FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL;
+
+ if ( mOrder < 10 )
+ {
+ mOrder = estimateBandPassOrder( mSampleRate, mPassFrequencyBegin, mPassFrequencyEnd,
+ mPassRipple, mStopRipple );
+ }
+
+ // Ensure even order since we're designing a Type 1 filter
+ if ( mOrder % 2 == 1 )
+ {
+ mOrder++;
+ }
+
+ FIRFilterSpecification spec = new FIRFilterSpecification( type, mOrder, mGridDensity );
+
+ //Set pass band to stop at about 80% of channel bandwidth as starting point
+ int passBandStop = (int)(mChannelBandwidth * .80);
+
+ //Set the stop band to start at 40% into the adjacent channel, since we're oversampling by 2.
+ int stopBandStart = (int)(mChannelBandwidth * 1.40);
+
+ FrequencyBand passBand = new FrequencyBand(mSampleRate, 0, passBandStop, 1.0, mPassRipple);
+ FrequencyBand stopBand = new FrequencyBand(mSampleRate, stopBandStart, (int)(mSampleRate / 2),
+ 0.0, mStopRipple);
+
+ spec.addFrequencyBand(passBand);
+ spec.addFrequencyBand(stopBand);
+
+ return spec;
+ }
+ }
+
+
+
+ // /**
// * Converts the set of band amplitude values into a set of weights
// * normalized to the largest band amplitude value
// */
@@ -918,7 +1074,7 @@ public class FIRFilterSpecification
*
* @return estimated filter length
*/
- private static int estimateFilterOrder( int sampleRate, int frequency1, int frequency2,
+ public static int estimateFilterOrder( int sampleRate, int frequency1, int frequency2,
double passBandRipple, double stopBandRipple )
{
double df = Math.abs( frequency2 - frequency1 ) / (double) sampleRate;
diff --git a/src/dsp/filter/fir/remez/PolyphaseChannelizerDesigner.java b/src/dsp/filter/fir/remez/PolyphaseChannelizerDesigner.java
new file mode 100644
index 00000000..f7436798
--- /dev/null
+++ b/src/dsp/filter/fir/remez/PolyphaseChannelizerDesigner.java
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * sdrtrunk
+ * Copyright (C) 2014-2017 Dennis Sheirer
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see
+ *
+ ******************************************************************************/
+package dsp.filter.fir.remez;
+
+import dsp.filter.fir.FIRFilterSpecification;
+
+public class PolyphaseChannelizerDesigner extends FIRFilterSpecification
+{
+ /**
+ * Constructs a Polyphase Channelizer Remez filter specification.
+ *
+ * Constructor is package private -- use the PolyphaseChannelizerFilterFactory to create a filter
+ *
+ * @param order of the filter
+ * @param gridDensity to use when calculating the coefficients
+ */
+ PolyphaseChannelizerDesigner(int order, int gridDensity)
+ {
+ super(FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL, order, gridDensity);
+ }
+}
diff --git a/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java b/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java
new file mode 100644
index 00000000..04417fa3
--- /dev/null
+++ b/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java
@@ -0,0 +1,132 @@
+/*******************************************************************************
+ * sdrtrunk
+ * Copyright (C) 2014-2017 Dennis Sheirer
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see
+ *
+ ******************************************************************************/
+package dsp.filter.fir.remez;
+
+import dsp.filter.fir.FIRFilterSpecification;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+public class PolyphaseChannelizerFilterFactory
+{
+ private final static Logger mLog = LoggerFactory.getLogger(PolyphaseChannelizerFilterFactory.class);
+
+ private static final double OBJECTIVE_BAND_EDGE_COEFFICIENT_AMPLITUDE = Math.sqrt(2.0) / 2.0; //.707xxx
+
+ public static float[] getFilter(int sampleRate, int channelBandwidth)
+ {
+ FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL;
+
+ int passBandStart = 0;
+ int gridDensity = 32;
+
+ //Set pass band to stop at about 80% of channel bandwidth as starting point
+// int passBandStop = (int)(channelBandwidth * 0.95);
+ int passBandStop = channelBandwidth;
+
+ //Set the stop band to start at 40% into the adjacent channel, since we're oversampling by 2.
+ int stopBandStart = (int)(channelBandwidth * 1.40);
+
+ double passRipple = 0.1d;
+ double stopRipple = 0.01d;
+
+// int order = FIRFilterSpecification.estimateBandPassOrder(sampleRate, passBandStart, passBandStop,
+// passRipple, stopRipple);
+ int order = FIRFilterSpecification.estimateFilterOrder( sampleRate, passBandStop,
+ stopBandStart, passRipple, stopRipple );
+
+ // Ensure even order since we're designing a Type 1 filter
+ if(order % 2 == 1)
+ {
+ order++;
+ }
+
+ FIRFilterSpecification specification = new PolyphaseChannelizerDesigner(order, gridDensity);
+
+ FIRFilterSpecification.FrequencyBand passBand = new FIRFilterSpecification.FrequencyBand(sampleRate, 0,
+ passBandStop, 1.0, passRipple);
+ FIRFilterSpecification.FrequencyBand stopBand = new FIRFilterSpecification.FrequencyBand(sampleRate,
+ stopBandStart, (int)(sampleRate / 2), 0.0, stopRipple);
+
+ specification.addFrequencyBand(passBand);
+ specification.addFrequencyBand(stopBand);
+
+ boolean complete = false;
+ int stepSize = 100;
+ double bandEdgeFrequency = Math.cos(Math.PI * (double)channelBandwidth / (double)(sampleRate / 2));
+ float[] filter = null;
+
+ while(!complete)
+ {
+ mLog.debug("Step Size: " + stepSize);
+
+ RemezFIRFilterDesigner2 designer = new RemezFIRFilterDesigner2(specification);
+
+ complete = true; //temporary
+
+ double bandEdgeAmplitude = designer.getFrequencyResponse(bandEdgeFrequency);
+ mLog.debug("Coefficient Amplitude at Band Edge is: " + bandEdgeAmplitude + " for frequency: " + passBandStop);
+
+ if(bandEdgeAmplitude < OBJECTIVE_BAND_EDGE_COEFFICIENT_AMPLITUDE)
+ {
+ if(stepSize <= 1)
+ {
+ complete = true;
+ try
+ {
+ filter = designer.getImpulseResponse();
+ }
+ catch(Exception e)
+ {
+ mLog.debug("Error getting filter response", e);
+ }
+ }
+ else
+ {
+ passBandStop += stepSize;
+ stepSize /= 10;
+ }
+ }
+ else
+ {
+ passBandStop -= stepSize;
+ }
+
+ if(!complete)
+ {
+ specification.clearFrequencyBands();
+
+ passBand = new FIRFilterSpecification.FrequencyBand(sampleRate, 0,
+ passBandStop, 1.0, passRipple);
+ specification.addFrequencyBand(passBand);
+ specification.addFrequencyBand(stopBand);
+ }
+ }
+
+ return filter;
+ }
+
+ public static void main(String[] args)
+ {
+ mLog.debug("Starting ...");
+
+ float[] filter = getFilter(2500000, 12500);
+
+ mLog.debug("Finished ...");
+ }
+}
diff --git a/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java
new file mode 100644
index 00000000..2890d939
--- /dev/null
+++ b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java
@@ -0,0 +1,707 @@
+/*******************************************************************************
+ * sdrtrunk
+ * Copyright (C) 2014-2017 Dennis Sheirer
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see
+ *
+ ******************************************************************************/
+package dsp.filter.fir.remez;
+
+import dsp.filter.design.FilterDesignException;
+import dsp.filter.fir.FIRFilterSpecification;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import java.math.BigDecimal;
+import java.math.RoundingMode;
+import java.util.ArrayList;
+import java.util.Iterator;
+import java.util.List;
+
+public class RemezFIRFilterDesigner2
+{
+ private final static Logger mLog = LoggerFactory.getLogger(RemezFIRFilterDesigner2.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 FIRFilterSpecification mSpecification;
+ private Grid mGrid;
+ private List mExtremalIndices;
+
+ private BigDecimal[] 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 RemezFIRFilterDesigner2(FIRFilterSpecification specification)
+ {
+ mSpecification = specification;
+
+ design();
+ }
+
+ /**
+ * Designs the filter according to the specification
+ */
+ private void design()
+ {
+ mLog.debug("Designing filter ...");
+ mGrid = new Grid(mSpecification);
+ mExtremalIndices = getInitialExtremalIndices();
+
+ int iterationCount = 0;
+
+ try
+ {
+ do
+ {
+ mLog.debug("Calculating grid frequency response");
+ calculateGridFrequencyResponse();
+ mLog.debug("Calculating grid error");
+ calculateGridError();
+ mLog.debug("Finding extremal indices");
+ findExtremalIndices();
+ mLog.debug("Checking Convergence");
+ checkConvergence();
+
+ iterationCount++;
+
+ mLog.debug("Iteration: " + 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
+ */
+ public double getFrequencyResponse(double cosineOfFrequency)
+ {
+ BigDecimal numerator = BigDecimal.ZERO;
+ BigDecimal denominator = BigDecimal.ZERO;
+
+ for(int k = 0; k < mExtremalIndices.size() - 1; k++)
+ {
+ BigDecimal cosineDelta = BigDecimal.valueOf(cosineOfFrequency -
+ mGrid.getCosineFrequencyGrid()[mExtremalIndices.get(k)]);
+
+ //If this frequency is close to one of the polynomial points, use c[k] for that point
+ if(cosineDelta.abs().compareTo(BigDecimal.valueOf(1.0e-7)) < 0)
+ {
+ numerator = BigDecimal.valueOf(mIdealFrequencyResponse[k]);
+ denominator = BigDecimal.ONE;
+ break;
+ }
+ else
+ {
+ BigDecimal dkOverCosineDelta = mD[k].divide(cosineDelta, 25, RoundingMode.HALF_UP);
+
+ numerator = numerator.add(dkOverCosineDelta.multiply(BigDecimal.valueOf(mIdealFrequencyResponse[k])));
+ denominator = denominator.add(dkOverCosineDelta);
+ }
+ }
+
+ return numerator.divide(denominator, 25, RoundingMode.HALF_UP).doubleValue();
+ }
+
+ /**
+ * 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 getInitialExtremalIndices()
+ {
+ int count = mSpecification.getExtremaCount();
+
+ int density = mSpecification.getGridDensity();
+
+ List 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()
+ {
+ mLog.debug("\tCalculating B");
+ BigDecimal[] b = calculateB();
+
+ mLog.debug("\tCalculating Delta");
+ calculateDelta(b);
+
+ mLog.debug("\tCalculating C");
+ calculateC();
+
+ mLog.debug("\tCalculating D");
+ calculateD(b);
+
+ mLog.debug("\tUpdating Grid Frequency Response");
+ 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++)
+ {
+ //TODO: parallelize this part!!!!
+ mGridFrequencyResponse[i] = getFrequencyResponse(gridFrequencyCosines[i]);
+
+ if(i % 1000 == 0)
+ {
+ mLog.debug("\t\tUpdating Frequency Response - index:" + i + " of " +
+ mGridFrequencyResponse.length + " Value:" + mGridFrequencyResponse[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 BigDecimal[] calculateB()
+ {
+ int length = mExtremalIndices.size();
+
+ BigDecimal[] b = new BigDecimal[length];
+
+ for(int k = 0; k < length; k++)
+ {
+ b[k] = BigDecimal.ONE;
+
+ 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] = b[k].multiply(BigDecimal.valueOf(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(BigDecimal[] b)
+ {
+ BigDecimal numerator = BigDecimal.ZERO;
+ BigDecimal denominator = BigDecimal.ZERO;
+
+ double sign = 1.0;
+
+ for(int k = 0; k < b.length; k++)
+ {
+ if(k < mExtremalIndices.size())
+ {
+ int extremalIndex = mExtremalIndices.get(k);
+
+ numerator = numerator.add(b[k].multiply(BigDecimal.valueOf(mGrid.getDesiredResponse()[extremalIndex])));
+ denominator = denominator.add(b[k].multiply(BigDecimal.valueOf(sign / mGrid.getWeight()[extremalIndex])));
+ sign = -sign;
+ }
+ else
+ {
+ mLog.error("Something went wrong -- the length of b exceeds the set of extremal indices");
+ }
+ }
+
+ mDelta = numerator.divide(denominator, 1023, BigDecimal.ROUND_HALF_DOWN).doubleValue();
+ }
+
+ /**
+ * 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(BigDecimal[] b)
+ {
+ int length = mExtremalIndices.size() - 1;
+
+ mD = new BigDecimal[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].multiply(BigDecimal.valueOf((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 indicesToRemove = new ArrayList<>();
+
+ Iterator 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;
+ }
+}