GitHub/sdrtrunk
1
0
Fork 0
mirror of https://github.com/DSheirer/sdrtrunk.git synced 2026-11-12 08:37:25 -05:00
Code Issues Packages Projects Releases 2 Wiki Activity

Adds a new remez designer that uses BigDecimal to handle the larger values inherent in longer length filters. However, this designer is hugely inefficient and takes forever to calculate the filter. More experimentation needed.

Browse Source
This commit is contained in:
Denny
2017-03-02 19:17:22 -05:00
parent cd4eaedfcc
commit 17c16e16ee
5 changed files with 1059 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
src/dsp/filter/design/FilterViewer.java
View File

@@ -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 )

172
src/dsp/filter/fir/FIRFilterSpecification.java
View File

@@ -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 <http://www.gnu.org/licenses/>
*
******************************************************************************/
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<FrequencyBand> 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;

37
src/dsp/filter/fir/remez/PolyphaseChannelizerDesigner.java Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>
*
******************************************************************************/
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);
}
}

132
src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>
*
******************************************************************************/
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 ...");
}
}

707
src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>
*
******************************************************************************/
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<Integer> 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<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()
{
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<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;
}
}