diff --git a/.idea/libraries/commons_math3_3_6_1.xml b/.idea/libraries/commons_math3_3_6_1.xml new file mode 100644 index 00000000..b692344c --- /dev/null +++ b/.idea/libraries/commons_math3_3_6_1.xml @@ -0,0 +1,9 @@ + + + + + + + + + \ No newline at end of file diff --git a/sdrtrunk.iml b/sdrtrunk.iml index 04cc0636..0621246d 100644 --- a/sdrtrunk.iml +++ b/sdrtrunk.iml @@ -21,5 +21,6 @@ + \ No newline at end of file diff --git a/src/dsp/filter/design/FilterViewer.java b/src/dsp/filter/design/FilterViewer.java index 348645b3..9bdf1ac0 100644 --- a/src/dsp/filter/design/FilterViewer.java +++ b/src/dsp/filter/design/FilterViewer.java @@ -2,7 +2,9 @@ package dsp.filter.design; import dsp.filter.FilterFactory; import dsp.filter.fir.FIRFilterSpecification; +import dsp.filter.fir.remez.PolyphaseChannelizerFilterFactory; import dsp.filter.fir.remez.RemezFIRFilterDesigner; +import dsp.filter.fir.remez.RemezFIRFilterDesignerWithLagrange; import javafx.application.Application; import javafx.scene.Scene; import javafx.stage.Stage; @@ -51,15 +53,15 @@ public class FilterViewer extends Application // .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 ) @@ -76,20 +78,22 @@ public class FilterViewer extends Application // 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); - } +// try +// { +// RemezFIRFilterDesignerWithLagrange designer = new RemezFIRFilterDesignerWithLagrange(specification); +//// RemezFIRFilterDesigner designer = new RemezFIRFilterDesigner(specification); +// +// if(designer.isValid()) +// { +// taps = designer.getImpulseResponse(); +// } +// } +// catch(FilterDesignException fde) +// { +// mLog.error("Filter design error", fde); +// } + taps = PolyphaseChannelizerFilterFactory.getFilter(10000000, 312500, 0.21); return taps; } } diff --git a/src/dsp/filter/fir/FIRFilterSpecification.java b/src/dsp/filter/fir/FIRFilterSpecification.java index ab599f89..da4aa40a 100644 --- a/src/dsp/filter/fir/FIRFilterSpecification.java +++ b/src/dsp/filter/fir/FIRFilterSpecification.java @@ -18,879 +18,677 @@ ******************************************************************************/ package dsp.filter.fir; +import dsp.filter.fir.remez.FIRLinearPhaseFilterType; + import java.text.DecimalFormat; import java.util.ArrayList; import java.util.Collection; import java.util.List; -import dsp.filter.fir.remez.FIRLinearPhaseFilterType; - public class FIRFilterSpecification { - private static final DecimalFormat DECIMAL_FORMATTER = new DecimalFormat( "0.00000" ); + private static final DecimalFormat DECIMAL_FORMATTER = new DecimalFormat("0.00000"); - private FIRLinearPhaseFilterType mRemezFilterType; - private int mOrder; - private int mGridDensity = 16; - private List mFrequencyBands = new ArrayList<>(); + private FIRLinearPhaseFilterType mRemezFilterType; + private int mOrder; + private int mGridDensity = 16; + private List mFrequencyBands = new ArrayList<>(); - /** - * 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 gridDensity to use when calculating the coefficients - */ - protected FIRFilterSpecification( FIRLinearPhaseFilterType type, int order, int gridDensity ) - { - mRemezFilterType = type; - mOrder = order; - mGridDensity = gridDensity; - } + /** + * 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 gridDensity to use when calculating the coefficients + */ + protected FIRFilterSpecification(FIRLinearPhaseFilterType type, int order, int gridDensity) + { + mRemezFilterType = type; + mOrder = order; + mGridDensity = gridDensity; + } - public String toString() - { - StringBuilder sb = new StringBuilder(); - sb.append( "Type: " ); - sb.append( mRemezFilterType.name() ); - sb.append( " Order: " ); - sb.append( mOrder ); - sb.append( " Extrema: " ); - sb.append( getExtremaCount() ); - sb.append( " Grid Density: " ); - sb.append( mGridDensity ); - sb.append( " Grid Size: " ); - sb.append( getGridSize() ); + public String toString() + { + StringBuilder sb = new StringBuilder(); + sb.append("Type: "); + sb.append(mRemezFilterType.name()); + sb.append(" Order: "); + sb.append(mOrder); + sb.append(" Extrema: "); + sb.append(getExtremaCount()); + sb.append(" Grid Density: "); + sb.append(mGridDensity); + sb.append(" Grid Size: "); + sb.append(getGridSize()); - for ( FrequencyBand band : mFrequencyBands ) - { - sb.append( "\n" ); - sb.append( band.toString() ); - } + for(FrequencyBand band : mFrequencyBands) + { + sb.append("\n"); + sb.append(band.toString()); + } - return sb.toString(); - } + return sb.toString(); + } - public void addFrequencyBand( FrequencyBand band ) - { - mFrequencyBands.add( band ); + public void addFrequencyBand(FrequencyBand band) + { + mFrequencyBands.add(band); - updateGridSize(); - } + updateGridSize(); + } - public void addFrequencyBands( Collection bands ) - { - mFrequencyBands.addAll( bands ); + public void addFrequencyBands(Collection bands) + { + mFrequencyBands.addAll(bands); - updateGridSize(); - } + updateGridSize(); + } /** * Removes all frequency bands from this filter specification */ - public void clearFrequencyBands() + public void clearFrequencyBands() { mFrequencyBands.clear(); } - public List getFrequencyBands() - { - return mFrequencyBands; - } - - /** - * Returns the total bandwidth of the defined frequency bands in the specification. - * - * @return total bandwidth -- should be less than or equal to 0.5 - */ - public double getTotalBandwidth() - { - double bandwidth = 0.0; - - for ( FrequencyBand band : mFrequencyBands ) - { - bandwidth += band.getBandWidth(); - } - - return bandwidth; - } - - /** - * Creates a Low-Pass FIR filter specification builder that allows you to define the filter - * parameters and create a filter specification. - */ - public static LowPassBuilder lowPassBuilder() - { - return new LowPassBuilder(); - } - - /** - * Creates a High-Pass FIR filter specification builder that allows you to define the filter - * parameters and create a filter specification. - */ - public static HighPassBuilder highPassBuilder() - { - return new HighPassBuilder(); - } - - /** - * Creates a Band-Pass FIR filter specification builder that allows you to define the filter - * parameters and create a filter specification. - */ - public static BandPassBuilder bandPassBuilder() - { - return new BandPassBuilder(); - } - - public FIRLinearPhaseFilterType getFilterType() - { - return mRemezFilterType; - } - - public boolean isSymmetrical() - { - return mRemezFilterType == FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL || - mRemezFilterType == FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL; - } - - public boolean isAntiSymmetrical() - { - return mRemezFilterType == FIRLinearPhaseFilterType.TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL || - mRemezFilterType == FIRLinearPhaseFilterType.TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL; - } - - public int getOrder() - { - return mOrder; - } - - public boolean isEvenOrder() - { - return mOrder % 2 == 0; - } - - public boolean isOddOrder() - { - return mOrder % 2 == 1; - } - - public int getFilterLength() - { - return mOrder + 1; - } - - public int getBandCount() - { - return mFrequencyBands.size(); - } - - /** - * Grid density setting. Default value is 16. - */ - public int getGridDensity() - { - return mGridDensity; - } - - /** - * Sets the grid density setting. Default value is 16. - */ - public void setGridDensity( int density ) - { - mGridDensity = density; - - updateGridSize(); - } - - /** - * Maximum number of alternations/extrema (L + 2 or M + 2) - */ - public int getExtremaCount() - { - return getHalfFilterOrder() + 2; - } - - /** - * Half of the filter order. - * - * This value is referred to as L or M, depending on the source. - */ - public int getHalfFilterOrder() - { - switch ( mRemezFilterType ) - { - case TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL: - return ( mOrder ) / 2; - case TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL: - return ( mOrder - 1 ) / 2; - case TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL: - return ( mOrder - 2 ) / 2; - case TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL: - return ( mOrder - 1 ) / 2; - default: - return 0; - } - } - - /** - * Determines the max band ripple from across the frequency bands in this specification. - */ - public double getMaxBandAmplitude() - { - double maxRippleAmplitude = 0.0; - - for ( FrequencyBand band : mFrequencyBands ) - { - double bandRippleAmplitude = band.getRippleAmplitude(); - - if ( bandRippleAmplitude > maxRippleAmplitude ) - { - maxRippleAmplitude = bandRippleAmplitude; - } - } - - return maxRippleAmplitude; - } - - public int getNominalGridSize() - { - return ( getExtremaCount() - 1 ) * getGridDensity() + 1; - } - - public int getGridSize() - { - int gridSize = 0; - - for ( FrequencyBand band : mFrequencyBands ) - { - gridSize += band.getGridSize(); - } - - return gridSize; - } - - private void updateGridSize() - { - if ( !mFrequencyBands.isEmpty() ) - { - int gridSize = getNominalGridSize(); - - double totalBandwidth = getTotalBandwidth(); - - for ( FrequencyBand band : mFrequencyBands ) - { - band.setGridSize( gridSize, totalBandwidth ); - } - } - } - - // /** - // * Determines the grid size needed to represent the percentage of weighted - // frequency bands - // * over the entire frequency band. - // */ - // public int getGridSize() - // { - // int gridSize = 0; - // - // int maxGridSize = getMaxGridSize(); - // - // //Add the proportional index counts from each weighted band - // for( FrequencyBand band: mFrequencyBands ) - // { - // gridSize += band.getBandGridSize( maxGridSize ); - // } - // - // if( getSymmetry() == Symmetry.NEGATIVE ) - // { - // gridSize--; - // } - // - // return gridSize; - // } - - /** - * Dense frequency spacing for extrema count * grid density over a frequency range of 0.0 - 1.0 - */ - public double getGridFrequencyInterval() - { - return getTotalBandwidth() / (double) ( ( getGridSize() - mFrequencyBands.size() ) ); - } - - // /** - // * ''' FIR order estimator (lowpass, highpass, bandpass, mulitiband). - // * - // * (n, fo, ao, w) = remezord (f, a, dev) (n, fo, ao, w) = remezord (f, a, - // * dev, fs) - // * - // * (n, fo, ao, w) = remezord (f, a, dev) finds the approximate order, - // * normalized frequency band edges, frequency band amplitudes, and weights - // * that meet input specifications f, a, and dev, to use with the remez - // * command. - // * - // * f is a sequence of frequency band edges (between 0 and Fs/2, where Fs - // is - // * the sampling frequency), and a is a sequence specifying the desired - // * amplitude on the bands defined by f. The length of f is twice the - // length - // * of a, minus 2. The desired function is piecewise constant. - // * - // * dev is a sequence the same size as a that specifies the maximum - // allowable - // * deviation or ripples between the frequency response and the desired - // * amplitude of the output filter, for each band. - // * - // * Use remez with the resulting order n, frequency sequence fo, amplitude - // * response sequence ao, and weights w to design the filter b which - // * approximately meets the specifications given by remezord input - // parameters - // * f, a, and dev: - // * - // * b = remez (n, fo, ao, w) - // * - // * (n, fo, ao, w) = remezord (f, a, dev, Fs) specifies a sampling - // frequency - // * Fs. - // * - // * Fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can - // * therefore specify band edges scaled to a particular applications - // sampling - // * frequency. - // * - // * In some cases remezord underestimates the order n. If the filter does - // not - // * meet the specifications, try a higher order such as n+1 or n+2. ''' - // * - // * @param fcuts - // * @param mags - // * @param devs - // * @param fsamp - // */ - // // TODO: rename to calculateFilterOrder() - // public static FIRFilterSpecification remezOrd( int[] fcutoffs, - // double[] mags, double[] devs, int fsamp ) - // { - // // def remezord (fcuts, mags, devs, fsamp = 2): - // - // double[] fcuts = new double[ fcutoffs.length ]; - // - // for ( int x = 0; x < fcuts.length; x++ ) - // { - // fcuts[ x ] = (double) fcutoffs[ x ] / (double) fsamp; - // } - // - // int nf = fcuts.length; - // int nm = mags.length; - // int nd = devs.length; - // int nbands = nm; - // - // if ( mags.length != devs.length ) - // { - // throw new IllegalArgumentException( - // "Length of magnitudes and deviations must be equal" ); - // } - // - // if ( fcuts.length != ( 2 * ( mags.length - 1 ) ) ) - // { - // throw new IllegalArgumentException( - // "Length of frequency band edges must be 2 * mags length - 2" ); - // } - // - // for ( int x = 0; x < mags.length; x++ ) - // { - // if ( mags[ x ] != 0 ) - // { - // devs[ x ] = devs[ x ] / mags[ x ]; - // } - // } - // - // // Separate the passband and stopband edges - // double[] f1 = new double[ fcuts.length / 2 ]; - // double[] f2 = new double[ fcuts.length / 2 ]; - // - // for ( int x = 0; x < fcuts.length / 2; x++ ) - // { - // f1[ x ] = fcuts[ 2 * x ]; - // f2[ x ] = fcuts[ 2 * x + 1 ]; - // } - // - // int n = 0; - // - // double min_delta = 2.0; - // - // for ( int x = 0; x < f1.length; x++ ) - // { - // if ( f2[ x ] - f1[ x ] < min_delta ) - // { - // n = x; - // min_delta = f2[ x ] - f1[ x ]; - // } - // } - // - // double l; - // - // if ( nbands == 2 ) - // { - // // lowpass or highpass case (use formula) - // l = lporder( f1[ n ], f2[ n ], devs[ 0 ], devs[ 1 ] ); - // } else - // { - // // bandpass or multipass case - // // try different lowpasses and take the worst one that - // // goes through the BP specs - // l = 0; - // - // for ( int x = 1; x < nbands - 1; x++ ) - // { - // double l1 = lporder( f1[ x - 1 ], f2[ x - 1 ], devs[ x ], - // devs[ x - 1 ] ); - // double l2 = lporder( f1[ x ], f2[ x ], devs[ x ], - // devs[ x + 1 ] ); - // l = Math.max( l1, l2 ); - // } - // } - // - // n = (int) ( Math.ceil( l ) - 1 ); // need order, not length for remez - // - // // // cook up remez compatible result - // double[] ff = new double[ fcuts.length + 2 ]; - // ff[ 0 ] = 0.0; - // ff[ ff.length - 1 ] = 1.0; - // - // System.arraycopy( fcuts, 0, ff, 1, fcuts.length ); - // - // for ( int x = 1; x < ff.length - 1; x++ ) - // { - // ff[ x ] *= 2; - // } - // - // double[] aa = new double[ mags.length * 2 ]; - // - // for ( int x = 0; x < mags.length; x++ ) - // { - // aa[ 2 * x ] = mags[ x ]; - // aa[ 2 * x + 1 ] = mags[ x ]; - // } - // - // double max_dev = 0; - // - // for ( double dev : devs ) - // { - // if ( dev > max_dev ) - // { - // max_dev = dev; - // } - // } - // - // double[] wts = new double[ devs.length ]; - // - // for ( int x = 0; x < devs.length; x++ ) - // { - // wts[ x ] = max_dev / devs[ x ]; - // } - // - // return null; - // // return new RemezFilterSpecification( RemezFilterType.BANDPASS, n, ff, - // // aa, wts, 16); - // } - - /** - * Low-pass FIR filter specification builder - */ - public static class LowPassBuilder - { - private int mSampleRate; - private int mOrder; - private int mGridDensity = 16; - private int mPassBandEndFrequency; - private int mStopBandStartFrequency; - private double mPassBandRipple; - private double mStopBandRipple; - private double mPassBandAmplitude = 1.0; - private double mStopBandAmplitude = 0.0; - - public LowPassBuilder() - { - } - - public LowPassBuilder sampleRate( int sampleRateHz ) - { - mSampleRate = sampleRateHz; - return this; - } - - /** - * Sets the filter order. If the order is not specified, it will be calculated from the - * other parameters. - */ - public LowPassBuilder order( int order ) - { - mOrder = order; - return this; - } - - public LowPassBuilder gridDensity( int density ) - { - mGridDensity = density; - return this; - } - - public LowPassBuilder passBandCutoff( int frequencyHz ) - { - mPassBandEndFrequency = frequencyHz; - return this; - } - - public LowPassBuilder stopBandStart( int frequencyHz ) - { - mStopBandStartFrequency = frequencyHz; - return this; - } - - public LowPassBuilder passBandRipple( double rippleDb ) - { - mPassBandRipple = rippleDb; - return this; - } - - public LowPassBuilder stopBandRipple( double rippleDb ) - { - mStopBandRipple = rippleDb; - return this; - } - - public LowPassBuilder passBandAmplitude( double amplitude ) - { - mPassBandAmplitude = amplitude; - return this; - } - - public LowPassBuilder stopBandAmplitude( double amplitude ) - { - mStopBandAmplitude = amplitude; - return this; - } - - public FIRFilterSpecification build() - { - if ( mOrder < 6 ) - { - mOrder = estimateFilterOrder( mSampleRate, mPassBandEndFrequency, - mStopBandStartFrequency, mPassBandRipple, mStopBandRipple ); - } - - FIRLinearPhaseFilterType type = ( mOrder % 2 == 0 - ? FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL - : FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL ); - - FIRFilterSpecification spec = new FIRFilterSpecification( type, mOrder, mGridDensity ); - - FrequencyBand passBand = new FrequencyBand( mSampleRate, 0, mPassBandEndFrequency, - mPassBandAmplitude, mPassBandRipple ); - spec.addFrequencyBand( passBand ); - - FrequencyBand stopBand = new FrequencyBand( mSampleRate, mStopBandStartFrequency, - mSampleRate / 2, mStopBandAmplitude, mStopBandRipple ); - spec.addFrequencyBand( stopBand ); - - return spec; - } - } - - /** - * Remez high-pass FIR filter specification builder - */ - public static class HighPassBuilder - { - private int mSampleRate; - private int mOrder = -1; - private int mGridDensity = 16; - private int mPassBandStartFrequency; - private int mStopBandEndFrequency; - private double mPassBandRipple; - private double mStopBandRipple; - private double mPassBandAmplitude = 1.0; - private double mStopBandAmplitude = 0.0; - - public HighPassBuilder() - { - } - - public HighPassBuilder sampleRate( int sampleRateHz ) - { - mSampleRate = sampleRateHz; - return this; - } - - /** - * Sets the filter order. If the order is not specified, it will be calculated from the - * other parameters. - */ - public HighPassBuilder order( int order ) - { - mOrder = order; - return this; - } - - public HighPassBuilder gridDensity( int density ) - { - mGridDensity = density; - return this; - } - - public HighPassBuilder passBandStart( int frequencyHz ) - { - mPassBandStartFrequency = frequencyHz; - return this; - } - - public HighPassBuilder stopBandCutoff( int frequencyHz ) - { - mStopBandEndFrequency = frequencyHz; - return this; - } - - public HighPassBuilder passBandRipple( double rippleDb ) - { - mPassBandRipple = rippleDb; - return this; - } - - public HighPassBuilder stopBandRipple( double rippleDb ) - { - mStopBandRipple = rippleDb; - return this; - } - - public HighPassBuilder passBandAmplitude( double amplitude ) - { - mPassBandAmplitude = amplitude; - return this; - } - - public HighPassBuilder stopBandAmplitude( double amplitude ) - { - mStopBandAmplitude = amplitude; - return this; - } - - public FIRFilterSpecification build() - { - if ( mOrder < 6 ) - { - mOrder = estimateFilterOrder( mSampleRate, mStopBandEndFrequency, - mPassBandStartFrequency, mPassBandRipple, mStopBandRipple ); - } - - // High pass is a Type 1 filter -- ensure order is even - if ( mOrder % 2 == 1 ) - { - mOrder++; - } - - FIRFilterSpecification spec = new FIRFilterSpecification( - FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL, mOrder, - mGridDensity ); - - FrequencyBand stopBand = new FrequencyBand( mSampleRate, 0, mStopBandEndFrequency, - mStopBandAmplitude, mStopBandRipple ); - spec.addFrequencyBand( stopBand ); - - FrequencyBand passBand = new FrequencyBand( mSampleRate, mPassBandStartFrequency, - mSampleRate / 2, mPassBandAmplitude, mPassBandRipple ); - spec.addFrequencyBand( passBand ); - - return spec; - } - } - - /** - * Band-pass or band-stop FIR filter specification builder - */ - public static class BandPassBuilder - { - private int mOrder; - private int mGridDensity = 16; - private int mSampleRate; - private int mStopFrequency1; - private int mPassFrequencyBegin; - private int mPassFrequencyEnd; - private int mStopFrequency2; - private double mStopRipple = 0.01; - private double mPassRipple = 0.01; - private double mStopAmplitude = 0.0; - private double mPassAmplitude = 1.0; - - public BandPassBuilder() - { - } - - /** - * 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 BandPassBuilder order( int order ) - { - mOrder = order; - return this; - } - - /** - * Sets the frequency grid density. Default is 16. - */ - public BandPassBuilder gridDensity( int density ) - { - mGridDensity = density; - return this; - } - - /** - * Sets the filter sample rate in hertz - */ - public BandPassBuilder sampleRate( int sampleRate ) - { - mSampleRate = sampleRate; - return this; - } - - /** - * Specifies the ending frequency of the first stop band - */ - public BandPassBuilder stopFrequency1( int stopFrequency ) - { - mStopFrequency1 = stopFrequency; - return this; - } - - /** - * Specifies the beginning frequency of the second stop band - */ - public BandPassBuilder stopFrequency2( int stopFrequency ) - { - mStopFrequency2 = stopFrequency; - return this; - } - - /** - * Specifies the pass band beginning frequency - */ - public BandPassBuilder passFrequencyBegin( int passFrequency ) - { - mPassFrequencyBegin = passFrequency; - return this; - } - - /** - * Specifies the pass band ending frequency - */ - public BandPassBuilder passFrequencyEnd( int passFrequency ) - { - mPassFrequencyEnd = passFrequency; - return this; - } - - /** - * Specifies the ripple in dB for both stop bands - */ - public BandPassBuilder stopRipple( double rippleDb ) - { - mStopRipple = rippleDb; - return this; - } - - /** - * Specifies the ripple in dB for the pass band - */ - public BandPassBuilder passRipple( double rippleDb ) - { - mPassRipple = rippleDb; - return this; - } - - /** - * Specifies the amplitude of the stop bands. Default is 0.0 - */ - public BandPassBuilder stopAmplitude( double amplitude ) - { - mStopAmplitude = amplitude; - return this; - } - - /** - * Specifies the amplitude of the pass band. Default is 1.0 - */ - public BandPassBuilder passAmplitude( double amplitude ) - { - mPassAmplitude = amplitude; - return this; - } - - /** - * Validates the frequency bands to test for overlap and completeness. - */ - private void validateFrequencyBands() - { - if ( mStopFrequency1 >= mPassFrequencyBegin ) - { - throw new IllegalArgumentException( "Stop band 1 ending frequency [" - + mStopFrequency1 + "] must be less than the pass band start frequency [" - + mPassFrequencyBegin + "]" ); - } - - if ( mPassFrequencyBegin >= mPassFrequencyEnd ) - { - throw new IllegalArgumentException( "Pass band begin frequency [" - + mPassFrequencyBegin + "] must be less than the pass band end frequency [" - + mPassFrequencyEnd + "]" ); - } - - if ( mPassFrequencyEnd >= mStopFrequency2 ) - { - throw new IllegalArgumentException( "Pass band end frequency [" + mPassFrequencyEnd - + "] must be less than stop band 2 beginning frequency [" + mStopFrequency2 - + "]" ); - } - - if ( mStopFrequency2 >= ( mSampleRate / 2 ) ) - { - throw new IllegalArgumentException( "Stop band 2 beginning frequency [" - + mStopFrequency2 + "] must be less than half of the sample rate [" - + ( mSampleRate / 2 ) + "]" ); - } - } - - public FIRFilterSpecification build() - { - validateFrequencyBands(); - - 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 ); - - spec.addFrequencyBand( new FrequencyBand( mSampleRate, 0, mStopFrequency1, - mStopAmplitude, mStopRipple ) ); - spec.addFrequencyBand( new FrequencyBand( mSampleRate, mPassFrequencyBegin, - mPassFrequencyEnd, mPassAmplitude, mPassRipple ) ); - spec.addFrequencyBand( new FrequencyBand( mSampleRate, mStopFrequency2, - (int) ( mSampleRate / 2 ), mStopAmplitude, mStopRipple ) ); - - return spec; - } - } + public List getFrequencyBands() + { + return mFrequencyBands; + } + + /** + * Returns the total bandwidth of the defined frequency bands in the specification. + * + * @return total bandwidth -- should be less than or equal to 0.5 + */ + public double getTotalBandwidth() + { + double bandwidth = 0.0; + + for(FrequencyBand band : mFrequencyBands) + { + bandwidth += band.getBandWidth(); + } + + return bandwidth; + } + + /** + * Creates a Low-Pass FIR filter specification builder that allows you to define the filter + * parameters and create a filter specification. + */ + public static LowPassBuilder lowPassBuilder() + { + return new LowPassBuilder(); + } + + /** + * Creates a High-Pass FIR filter specification builder that allows you to define the filter + * parameters and create a filter specification. + */ + public static HighPassBuilder highPassBuilder() + { + return new HighPassBuilder(); + } + + /** + * Creates a Band-Pass FIR filter specification builder that allows you to define the filter + * parameters and create a filter specification. + */ + public static BandPassBuilder bandPassBuilder() + { + return new BandPassBuilder(); + } + + public FIRLinearPhaseFilterType getFilterType() + { + return mRemezFilterType; + } + + public boolean isSymmetrical() + { + return mRemezFilterType == FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL || + mRemezFilterType == FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL; + } + + public boolean isAntiSymmetrical() + { + return mRemezFilterType == FIRLinearPhaseFilterType.TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL || + mRemezFilterType == FIRLinearPhaseFilterType.TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL; + } + + public int getOrder() + { + return mOrder; + } + + public boolean isEvenOrder() + { + return mOrder % 2 == 0; + } + + public boolean isOddOrder() + { + return mOrder % 2 == 1; + } + + public int getFilterLength() + { + return mOrder + 1; + } + + public int getBandCount() + { + return mFrequencyBands.size(); + } + + /** + * Grid density setting. Default value is 16. + */ + public int getGridDensity() + { + return mGridDensity; + } + + /** + * Sets the grid density setting. Default value is 16. + */ + public void setGridDensity(int density) + { + mGridDensity = density; + + updateGridSize(); + } + + /** + * Maximum number of alternations/extrema (L + 2 or M + 2) + */ + public int getExtremaCount() + { + return getHalfFilterOrder() + 2; + } + + /** + * Half of the filter order. + * + * This value is referred to as L or M, depending on the source. + */ + public int getHalfFilterOrder() + { + switch(mRemezFilterType) + { + case TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL: + return (mOrder) / 2; + case TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL: + return (mOrder - 1) / 2; + case TYPE_3_ODD_LENGTH_EVEN_ORDER_ANTI_SYMMETRICAL: + return (mOrder - 2) / 2; + case TYPE_4_EVEN_LENGTH_ODD_ORDER_ANTI_SYMMETRICAL: + return (mOrder - 1) / 2; + default: + return 0; + } + } + + /** + * Determines the max band ripple from across the frequency bands in this specification. + */ + public double getMaxBandAmplitude() + { + double maxRippleAmplitude = 0.0; + + for(FrequencyBand band : mFrequencyBands) + { + double bandRippleAmplitude = band.getRippleAmplitude(); + + if(bandRippleAmplitude > maxRippleAmplitude) + { + maxRippleAmplitude = bandRippleAmplitude; + } + } + + return maxRippleAmplitude; + } + + public int getNominalGridSize() + { + return (getExtremaCount() - 1) * getGridDensity() + 1; + } + + public int getGridSize() + { + int gridSize = 0; + + for(FrequencyBand band : mFrequencyBands) + { + gridSize += band.getGridSize(); + } + + return gridSize; + } + + private void updateGridSize() + { + if(!mFrequencyBands.isEmpty()) + { + int gridSize = getNominalGridSize(); + + double totalBandwidth = getTotalBandwidth(); + + for(FrequencyBand band : mFrequencyBands) + { + band.setGridSize(gridSize, totalBandwidth); + } + } + } + + /** + * Dense frequency spacing for extrema count * grid density over a frequency range of 0.0 - 1.0 + */ + public double getGridFrequencyInterval() + { + return getTotalBandwidth() / (double)((getGridSize() - mFrequencyBands.size())); + } + + /** + * Low-pass FIR filter specification builder + */ + public static class LowPassBuilder + { + private int mSampleRate; + private int mOrder; + private int mGridDensity = 16; + private int mPassBandEndFrequency; + private int mStopBandStartFrequency; + private double mPassBandRipple; + private double mStopBandRipple; + private double mPassBandAmplitude = 1.0; + private double mStopBandAmplitude = 0.0; + + public LowPassBuilder() + { + } + + public LowPassBuilder sampleRate(int sampleRateHz) + { + mSampleRate = sampleRateHz; + return this; + } + + /** + * Sets the filter order. If the order is not specified, it will be calculated from the + * other parameters. + */ + public LowPassBuilder order(int order) + { + mOrder = order; + return this; + } + + public LowPassBuilder gridDensity(int density) + { + mGridDensity = density; + return this; + } + + public LowPassBuilder passBandCutoff(int frequencyHz) + { + mPassBandEndFrequency = frequencyHz; + return this; + } + + public LowPassBuilder stopBandStart(int frequencyHz) + { + mStopBandStartFrequency = frequencyHz; + return this; + } + + public LowPassBuilder passBandRipple(double rippleDb) + { + mPassBandRipple = rippleDb; + return this; + } + + public LowPassBuilder stopBandRipple(double rippleDb) + { + mStopBandRipple = rippleDb; + return this; + } + + public LowPassBuilder passBandAmplitude(double amplitude) + { + mPassBandAmplitude = amplitude; + return this; + } + + public LowPassBuilder stopBandAmplitude(double amplitude) + { + mStopBandAmplitude = amplitude; + return this; + } + + public FIRFilterSpecification build() + { + if(mOrder < 6) + { + mOrder = estimateFilterOrder(mSampleRate, mPassBandEndFrequency, + mStopBandStartFrequency, mPassBandRipple, mStopBandRipple); + } + + FIRLinearPhaseFilterType type = (mOrder % 2 == 0 + ? FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL + : FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL); + + FIRFilterSpecification spec = new FIRFilterSpecification(type, mOrder, mGridDensity); + + FrequencyBand passBand = new FrequencyBand(mSampleRate, 0, mPassBandEndFrequency, + mPassBandAmplitude, mPassBandRipple); + spec.addFrequencyBand(passBand); + + FrequencyBand stopBand = new FrequencyBand(mSampleRate, mStopBandStartFrequency, + mSampleRate / 2, mStopBandAmplitude, mStopBandRipple); + spec.addFrequencyBand(stopBand); + + return spec; + } + } + + /** + * Remez high-pass FIR filter specification builder + */ + public static class HighPassBuilder + { + private int mSampleRate; + private int mOrder = -1; + private int mGridDensity = 16; + private int mPassBandStartFrequency; + private int mStopBandEndFrequency; + private double mPassBandRipple; + private double mStopBandRipple; + private double mPassBandAmplitude = 1.0; + private double mStopBandAmplitude = 0.0; + + public HighPassBuilder() + { + } + + public HighPassBuilder sampleRate(int sampleRateHz) + { + mSampleRate = sampleRateHz; + return this; + } + + /** + * Sets the filter order. If the order is not specified, it will be calculated from the + * other parameters. + */ + public HighPassBuilder order(int order) + { + mOrder = order; + return this; + } + + public HighPassBuilder gridDensity(int density) + { + mGridDensity = density; + return this; + } + + public HighPassBuilder passBandStart(int frequencyHz) + { + mPassBandStartFrequency = frequencyHz; + return this; + } + + public HighPassBuilder stopBandCutoff(int frequencyHz) + { + mStopBandEndFrequency = frequencyHz; + return this; + } + + public HighPassBuilder passBandRipple(double rippleDb) + { + mPassBandRipple = rippleDb; + return this; + } + + public HighPassBuilder stopBandRipple(double rippleDb) + { + mStopBandRipple = rippleDb; + return this; + } + + public HighPassBuilder passBandAmplitude(double amplitude) + { + mPassBandAmplitude = amplitude; + return this; + } + + public HighPassBuilder stopBandAmplitude(double amplitude) + { + mStopBandAmplitude = amplitude; + return this; + } + + public FIRFilterSpecification build() + { + if(mOrder < 6) + { + mOrder = estimateFilterOrder(mSampleRate, mStopBandEndFrequency, + mPassBandStartFrequency, mPassBandRipple, mStopBandRipple); + } + + // High pass is a Type 1 filter -- ensure order is even + if(mOrder % 2 == 1) + { + mOrder++; + } + + FIRFilterSpecification spec = new FIRFilterSpecification( + FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL, mOrder, + mGridDensity); + + FrequencyBand stopBand = new FrequencyBand(mSampleRate, 0, mStopBandEndFrequency, + mStopBandAmplitude, mStopBandRipple); + spec.addFrequencyBand(stopBand); + + FrequencyBand passBand = new FrequencyBand(mSampleRate, mPassBandStartFrequency, + mSampleRate / 2, mPassBandAmplitude, mPassBandRipple); + spec.addFrequencyBand(passBand); + + return spec; + } + } + + /** + * Band-pass or band-stop FIR filter specification builder + */ + public static class BandPassBuilder + { + private int mOrder; + private int mGridDensity = 16; + private int mSampleRate; + private int mStopFrequency1; + private int mPassFrequencyBegin; + private int mPassFrequencyEnd; + private int mStopFrequency2; + private double mStopRipple = 0.01; + private double mPassRipple = 0.01; + private double mStopAmplitude = 0.0; + private double mPassAmplitude = 1.0; + + public BandPassBuilder() + { + } + + /** + * 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 BandPassBuilder order(int order) + { + mOrder = order; + return this; + } + + /** + * Sets the frequency grid density. Default is 16. + */ + public BandPassBuilder gridDensity(int density) + { + mGridDensity = density; + return this; + } + + /** + * Sets the filter sample rate in hertz + */ + public BandPassBuilder sampleRate(int sampleRate) + { + mSampleRate = sampleRate; + return this; + } + + /** + * Specifies the ending frequency of the first stop band + */ + public BandPassBuilder stopFrequency1(int stopFrequency) + { + mStopFrequency1 = stopFrequency; + return this; + } + + /** + * Specifies the beginning frequency of the second stop band + */ + public BandPassBuilder stopFrequency2(int stopFrequency) + { + mStopFrequency2 = stopFrequency; + return this; + } + + /** + * Specifies the pass band beginning frequency + */ + public BandPassBuilder passFrequencyBegin(int passFrequency) + { + mPassFrequencyBegin = passFrequency; + return this; + } + + /** + * Specifies the pass band ending frequency + */ + public BandPassBuilder passFrequencyEnd(int passFrequency) + { + mPassFrequencyEnd = passFrequency; + return this; + } + + /** + * Specifies the ripple in dB for both stop bands + */ + public BandPassBuilder stopRipple(double rippleDb) + { + mStopRipple = rippleDb; + return this; + } + + /** + * Specifies the ripple in dB for the pass band + */ + public BandPassBuilder passRipple(double rippleDb) + { + mPassRipple = rippleDb; + return this; + } + + /** + * Specifies the amplitude of the stop bands. Default is 0.0 + */ + public BandPassBuilder stopAmplitude(double amplitude) + { + mStopAmplitude = amplitude; + return this; + } + + /** + * Specifies the amplitude of the pass band. Default is 1.0 + */ + public BandPassBuilder passAmplitude(double amplitude) + { + mPassAmplitude = amplitude; + return this; + } + + /** + * Validates the frequency bands to test for overlap and completeness. + */ + private void validateFrequencyBands() + { + if(mStopFrequency1 >= mPassFrequencyBegin) + { + throw new IllegalArgumentException("Stop band 1 ending frequency [" + + mStopFrequency1 + "] must be less than the pass band start frequency [" + + mPassFrequencyBegin + "]"); + } + + if(mPassFrequencyBegin >= mPassFrequencyEnd) + { + throw new IllegalArgumentException("Pass band begin frequency [" + + mPassFrequencyBegin + "] must be less than the pass band end frequency [" + + mPassFrequencyEnd + "]"); + } + + if(mPassFrequencyEnd >= mStopFrequency2) + { + throw new IllegalArgumentException("Pass band end frequency [" + mPassFrequencyEnd + + "] must be less than stop band 2 beginning frequency [" + mStopFrequency2 + + "]"); + } + + if(mStopFrequency2 >= (mSampleRate / 2)) + { + throw new IllegalArgumentException("Stop band 2 beginning frequency [" + + mStopFrequency2 + "] must be less than half of the sample rate [" + + (mSampleRate / 2) + "]"); + } + } + + public FIRFilterSpecification build() + { + validateFrequencyBands(); + + 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); + + spec.addFrequencyBand(new FrequencyBand(mSampleRate, 0, mStopFrequency1, + mStopAmplitude, mStopRipple)); + spec.addFrequencyBand(new FrequencyBand(mSampleRate, mPassFrequencyBegin, + mPassFrequencyEnd, mPassAmplitude, mPassRipple)); + spec.addFrequencyBand(new FrequencyBand(mSampleRate, mStopFrequency2, + (int)(mSampleRate / 2), mStopAmplitude, mStopRipple)); + + return spec; + } + } /** * Builder for creating polyphase channelizer low-pass prototype filter. @@ -916,7 +714,7 @@ public class FIRFilterSpecification /** * Sets the filter sample rate in hertz */ - public PolyphaseChannelizerBuilder sampleRate(int sampleRate ) + public PolyphaseChannelizerBuilder sampleRate(int sampleRate) { mSampleRate = sampleRate; return this; @@ -936,7 +734,7 @@ public class FIRFilterSpecification * 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 ) + public PolyphaseChannelizerBuilder order(int order) { mOrder = order; return this; @@ -945,7 +743,7 @@ public class FIRFilterSpecification /** * Sets the frequency grid density. Default is 16. */ - public PolyphaseChannelizerBuilder gridDensity(int density ) + public PolyphaseChannelizerBuilder gridDensity(int density) { mGridDensity = density; return this; @@ -954,7 +752,7 @@ public class FIRFilterSpecification /** * Specifies the ripple in dB for both stop bands */ - public PolyphaseChannelizerBuilder stopRipple(double rippleDb ) + public PolyphaseChannelizerBuilder stopRipple(double rippleDb) { mStopRipple = rippleDb; return this; @@ -963,7 +761,7 @@ public class FIRFilterSpecification /** * Specifies the ripple in dB for the pass band */ - public PolyphaseChannelizerBuilder passRipple(double rippleDb ) + public PolyphaseChannelizerBuilder passRipple(double rippleDb) { mPassRipple = rippleDb; return this; @@ -972,7 +770,7 @@ public class FIRFilterSpecification /** * Specifies the amplitude of the stop bands. Default is 0.0 */ - public PolyphaseChannelizerBuilder stopAmplitude(double amplitude ) + public PolyphaseChannelizerBuilder stopAmplitude(double amplitude) { mStopAmplitude = amplitude; return this; @@ -981,7 +779,7 @@ public class FIRFilterSpecification /** * Specifies the amplitude of the pass band. Default is 1.0 */ - public PolyphaseChannelizerBuilder passAmplitude(double amplitude ) + public PolyphaseChannelizerBuilder passAmplitude(double amplitude) { mPassAmplitude = amplitude; return this; @@ -991,19 +789,19 @@ public class FIRFilterSpecification { FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL; - if ( mOrder < 10 ) + if(mOrder < 10) { - mOrder = estimateBandPassOrder( mSampleRate, mPassFrequencyBegin, mPassFrequencyEnd, - mPassRipple, mStopRipple ); + mOrder = estimateBandPassOrder(mSampleRate, mPassFrequencyBegin, mPassFrequencyEnd, + mPassRipple, mStopRipple); } // Ensure even order since we're designing a Type 1 filter - if ( mOrder % 2 == 1 ) + if(mOrder % 2 == 1) { mOrder++; } - FIRFilterSpecification spec = new FIRFilterSpecification( type, mOrder, mGridDensity ); + 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); @@ -1023,265 +821,273 @@ public class FIRFilterSpecification } - // /** - // * Converts the set of band amplitude values into a set of weights - // * normalized to the largest band amplitude value - // */ - // private static double[] getBandWeights( double... amplitudes ) - // { - // double maximum = 0.0; - // - // for ( double amplitude : amplitudes ) - // { - // if ( amplitude > maximum ) - // { - // maximum = amplitude; - // } - // } - // - // double[] weights = new double[ amplitudes.length ]; - // - // for ( int x = 0; x < amplitudes.length; x++ ) - // { - // weights[ x ] = maximum / amplitudes[ x ]; - // } - // - // return weights; - // } + // * Converts the set of band amplitude values into a set of weights + // * normalized to the largest band amplitude value + // */ + // private static double[] getBandWeights( double... amplitudes ) + // { + // double maximum = 0.0; + // + // for ( double amplitude : amplitudes ) + // { + // if ( amplitude > maximum ) + // { + // maximum = amplitude; + // } + // } + // + // double[] weights = new double[ amplitudes.length ]; + // + // for ( int x = 0; x < amplitudes.length; x++ ) + // { + // weights[ x ] = maximum / amplitudes[ x ]; + // } + // + // return weights; + // } - /** - * Calculates filter order from filter length - */ - private static int getFilterOrder( double length ) - { - return (int) ( Math.ceil( length ) - 1 ); - } + /** + * Calculates filter order from filter length + */ + private static int getFilterOrder(double length) + { + return (int)(Math.ceil(length) - 1); + } - /** - * Estimates filter length for a two-band (high or low pass) FIR filter. - * - * Note: this estimator doesn't work well if the transition frequency is near 0 or sampleRate/2 - * - * From Herrmann et al (1973), Practical design rules for optimum finite impulse response - * filters. Bell System Technical J., 52, 769-99 - * - * @param sampleRate in hertz - * @param frequency1 in hertz - * @param frequency2 in hertz - * @param passBandRipple pass band ripple in dB - * @param stopBandRipple stop band ripple in dB - * - * @return estimated filter length - */ - public static int estimateFilterOrder( int sampleRate, int frequency1, int frequency2, - double passBandRipple, double stopBandRipple ) - { - double df = Math.abs( frequency2 - frequency1 ) / (double) sampleRate; + /** + * Estimates filter length for a two-band (high or low pass) FIR filter. + * + * Note: this estimator doesn't work well if the transition frequency is near 0 or sampleRate/2 + * + * From Herrmann et al (1973), Practical design rules for optimum finite impulse response + * filters. Bell System Technical J., 52, 769-99 + * + * @param sampleRate in hertz + * @param frequency1 in hertz + * @param frequency2 in hertz + * @param passBandRipple pass band ripple in dB + * @param stopBandRipple stop band ripple in dB + * @return estimated filter length + */ + public static int estimateFilterOrder(int sampleRate, int frequency1, int frequency2, + double passBandRipple, double stopBandRipple) + { + double df = Math.abs(frequency2 - frequency1) / (double)sampleRate; - double ddp = Math - .log10( stopBandRipple <= passBandRipple ? passBandRipple : stopBandRipple ); - double dds = Math - .log10( stopBandRipple <= passBandRipple ? stopBandRipple : passBandRipple ); + double ddp = Math + .log10(stopBandRipple <= passBandRipple ? passBandRipple : stopBandRipple); + double dds = Math + .log10(stopBandRipple <= passBandRipple ? stopBandRipple : passBandRipple); - double a1 = 5.309e-3; - double a2 = 7.114e-2; - double a3 = -4.761e-1; - double a4 = -2.66e-3; - double a5 = -5.941e-1; - double a6 = -4.278e-1; + double a1 = 5.309e-3; + double a2 = 7.114e-2; + double a3 = -4.761e-1; + double a4 = -2.66e-3; + double a5 = -5.941e-1; + double a6 = -4.278e-1; - double b1 = 11.01217; - double b2 = 0.5124401; + double b1 = 11.01217; + double b2 = 0.5124401; - double t1 = a1 * ddp * ddp; - double t2 = a2 * ddp; - double t3 = a4 * ddp * ddp; - double t4 = a5 * ddp; + double t1 = a1 * ddp * ddp; + double t2 = a2 * ddp; + double t3 = a4 * ddp * ddp; + double t4 = a5 * ddp; - double dinf = ( ( t1 + t2 + a3 ) * dds ) + ( t3 + t4 + a6 ); - double ff = b1 + b2 * ( ddp - dds ); - double n = dinf / df - ff * df + 1.0; + double dinf = ((t1 + t2 + a3) * dds) + (t3 + t4 + a6); + double ff = b1 + b2 * (ddp - dds); + double n = dinf / df - ff * df + 1.0; - return (int) Math.ceil( n ); - } + return (int)Math.ceil(n); + } - /** - * Estimates band pass filter length; - * - * @param passBandStart frequency normalized to sample rate - * @param passBandEnd frequency normalized to sample rate - * @param passBandRippleDb - * @param stopBandRippleDb - */ - public static int estimateBandPassOrder( int sampleRate, int passBandStart, int passBandEnd, - double passBandRippleDb, double stopBandRippleDb ) - { - double df = (double) Math.abs( passBandEnd - passBandStart ) / (double) sampleRate; - double ddp = (double) Math.log10( passBandRippleDb ); - double dds = (double) Math.log10( stopBandRippleDb ); + /** + * Estimates band pass filter length; + * + * @param passBandStart frequency normalized to sample rate + * @param passBandEnd frequency normalized to sample rate + * @param passBandRippleDb + * @param stopBandRippleDb + */ + public static int estimateBandPassOrder(int sampleRate, int passBandStart, int passBandEnd, + double passBandRippleDb, double stopBandRippleDb) + { + double df = (double)Math.abs(passBandEnd - passBandStart) / (double)sampleRate; + double ddp = (double)Math.log10(passBandRippleDb); + double dds = (double)Math.log10(stopBandRippleDb); - double a1 = 0.01201; - double a2 = 0.09664; - double a3 = -0.51325; - double a4 = 0.00203; - double a5 = -0.57054; - double a6 = -0.44314; + double a1 = 0.01201; + double a2 = 0.09664; + double a3 = -0.51325; + double a4 = 0.00203; + double a5 = -0.57054; + double a6 = -0.44314; - double t1 = a1 * ddp * ddp; - double t2 = a2 * ddp; - double t3 = a4 * ddp * ddp; - double t4 = a5 * ddp; + double t1 = a1 * ddp * ddp; + double t2 = a2 * ddp; + double t3 = a4 * ddp * ddp; + double t4 = a5 * ddp; - double cinf = dds * ( t1 + t2 + a3 ) + t3 + t4 + a6; - double ginf = -14.6f * (double) Math.log10( passBandRippleDb / stopBandRippleDb ) - 16.9; - double n = cinf / df + ginf * df + 1.0; + double cinf = dds * (t1 + t2 + a3) + t3 + t4 + a6; + double ginf = -14.6f * (double)Math.log10(passBandRippleDb / stopBandRippleDb) - 16.9; + double n = cinf / df + ginf * df + 1.0; - return (int) Math.ceil( n ); - } + return (int)Math.ceil(n); + } - /** - * Frequency band - */ - public static class FrequencyBand - { - private int mGridSize; - private double mStart; - private double mEnd; - private double mAmplitude; - private double mRippleDB; + /** + * Frequency band + */ + public static class FrequencyBand + { + private int mGridSize; + private double mStart; + private double mEnd; + private double mAmplitude; + private double mRippleDB; + private Double mWeight; - /** - * Constructs a frequency band where the start and end frequencies are normalized to a - * sample rate of 1. - * - * @param start edge of the frequency band (0.0 - 1.0) - * @param end edge of the frequency band (0.0 - 1.0) - * @param amplitude of the frequency band on the unity scale (0.0 - 1.0) - * @param ripple of the frequency band specified as db ripple - */ - public FrequencyBand( double start, double end, double amplitude, double ripple ) - { - assert ( 0.0 <= mStart ); - assert ( mStart < mEnd ); - assert ( mEnd <= 1.0 ); - assert ( 0.0 <= amplitude && amplitude <= 1.0 ); - assert ( 0.0 <= ripple ); + /** + * Constructs a frequency band where the start and end frequencies are normalized to a + * sample rate of 1. + * + * @param start edge of the frequency band (0.0 - 1.0) + * @param end edge of the frequency band (0.0 - 1.0) + * @param amplitude of the frequency band on the unity scale (0.0 - 1.0) + * @param ripple of the frequency band specified as db ripple + */ + public FrequencyBand(double start, double end, double amplitude, double ripple) + { + assert (0.0 <= mStart); + assert (mStart <= mEnd); + assert (mEnd <= 1.0); + assert (0.0 <= amplitude && amplitude <= 1.0); + assert (0.0 <= ripple); - mStart = start; - mEnd = end; - mAmplitude = amplitude; - mRippleDB = ripple; - } + mStart = start; + mEnd = end; + mAmplitude = amplitude; + mRippleDB = ripple; + } - /** - * Constructs a frequency band for the sample rate and start and end frequencies. - * - * @param sampleRate in hertz - * @param start edge of the frequency band in hertz - * @param end edge of the frequency band in hertz - * @param amplitude of the frequency band relative to unity (0.0 - 1.0) - * @param ripple of the frequency band specified as db ripple - */ - public FrequencyBand( int sampleRate, int start, int end, double amplitude, double ripple ) - { - this( (double) start / (double) sampleRate, (double) end / (double) sampleRate, - amplitude, ripple ); - } + /** + * Constructs a frequency band for the sample rate and start and end frequencies. + * + * @param sampleRate in hertz + * @param start edge of the frequency band in hertz + * @param end edge of the frequency band in hertz + * @param amplitude of the frequency band relative to unity (0.0 - 1.0) + * @param ripple of the frequency band specified as db ripple + */ + public FrequencyBand(int sampleRate, int start, int end, double amplitude, double ripple) + { + this((double)start / (double)sampleRate, (double)end / (double)sampleRate, + amplitude, ripple); + } - public int getGridSize() - { - return mGridSize; - } + public FrequencyBand(int sampleRate, int start, int end, double amplitude, double ripple, double weight) + { + this(sampleRate, start, end, amplitude, ripple); + mWeight = weight; + } - public void setGridSize( int totalGridSize, double totalBandwidth ) - { - mGridSize = (int) Math - .ceil( (double) totalGridSize * ( getBandWidth() / totalBandwidth ) ); - } + public int getGridSize() + { + return mGridSize; + } - public String toString() - { - StringBuilder sb = new StringBuilder(); + public void setGridSize(int totalGridSize, double totalBandwidth) + { + mGridSize = Math.max(1, (int)Math.ceil((double)totalGridSize * (getBandWidth() / totalBandwidth))); + } - sb.append( "Band Start: " ); - sb.append( mStart ); - sb.append( " End:" ); - sb.append( mEnd ); - sb.append( " Amplitude: " ); - sb.append( mAmplitude ); - sb.append( " Ripple dB: " ); - sb.append( mRippleDB ); - sb.append( " Ripple Amplitude: " ); - sb.append( DECIMAL_FORMATTER.format( getRippleAmplitude() ) ); - sb.append( " Grid Size: " ); - sb.append( mGridSize ); + public String toString() + { + StringBuilder sb = new StringBuilder(); - return sb.toString(); - } + sb.append("Band Start: "); + sb.append(mStart); + sb.append(" End:"); + sb.append(mEnd); + sb.append(" Amplitude: "); + sb.append(mAmplitude); + sb.append(" Ripple dB: "); + sb.append(mRippleDB); + sb.append(" Ripple Amplitude: "); + sb.append(DECIMAL_FORMATTER.format(getRippleAmplitude())); + sb.append(" Grid Size: "); + sb.append(mGridSize); - /** - * Start frequency edge for this band normalized to 1 Hz. - */ - public double getStart() - { - return mStart; - } + return sb.toString(); + } - /** - * End frequency edge for this band normalized to 1 Hz. - */ - public double getEnd() - { - return mEnd; - } + /** + * Start frequency edge for this band normalized to 1 Hz. + */ + public double getStart() + { + return mStart; + } - /** - * Desired amplitude response for this frequency band - * - * @return amplitude (0.0 - 1.0) - */ - public double getAmplitude() - { - return mAmplitude; - } + /** + * End frequency edge for this band normalized to 1 Hz. + */ + public double getEnd() + { + return mEnd; + } - /** - * Frequency band ripple in decibels - */ - public double getRippleDB() - { - return mRippleDB; - } + /** + * Desired amplitude response for this frequency band + * + * @return amplitude (0.0 - 1.0) + */ + public double getAmplitude() + { + return mAmplitude; + } - /** - * Frequency band ripple as an amplitude value - */ - public double getRippleAmplitude() - { - return ( Math.pow( 10.0, ( mRippleDB / 20 ) ) - 1 ) - / ( Math.pow( 10.0, ( mRippleDB / 20 ) ) + 1 ); - } + /** + * Frequency band ripple in decibels + */ + public double getRippleDB() + { + return mRippleDB; + } - /** - * Calculates the weighting of this frequency band's ripple relative to the max ripple for a - * filter specification. - * - * @param maxRippleAmplitude for the filter specification - * - * @return weight of this frequency band - */ - public double getWeight( double maxRippleAmplitude ) - { - return 1.0 / ( getRippleAmplitude() / maxRippleAmplitude ); - } + /** + * Frequency band ripple as an amplitude value + */ + public double getRippleAmplitude() + { + return (Math.pow(10.0, (mRippleDB / 20)) - 1) + / (Math.pow(10.0, (mRippleDB / 20)) + 1); + } - /** - * Double-sided bandwidth of this frequency band - */ - public double getBandWidth() - { - return mEnd - mStart; - } - } + /** + * Calculates the weighting of this frequency band's ripple relative to the max ripple for a + * filter specification. + * + * @param maxRippleAmplitude for the filter specification + * @return weight of this frequency band + */ + public double getWeight(double maxRippleAmplitude) + { + if(mWeight == null) + { + return 1.0 / (getRippleAmplitude() / maxRippleAmplitude); + } + + return mWeight; + } + + /** + * Double-sided bandwidth of this frequency band + */ + public double getBandWidth() + { + return mEnd - mStart; + } + } } diff --git a/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java b/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java index 04417fa3..7e21f6fa 100644 --- a/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java +++ b/src/dsp/filter/fir/remez/PolyphaseChannelizerFilterFactory.java @@ -28,27 +28,25 @@ public class PolyphaseChannelizerFilterFactory private static final double OBJECTIVE_BAND_EDGE_COEFFICIENT_AMPLITUDE = Math.sqrt(2.0) / 2.0; //.707xxx - public static float[] getFilter(int sampleRate, int channelBandwidth) + public static float[] getFilter(int sampleRate, int channelBandwidth, double alpha) { FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL; int passBandStart = 0; - int gridDensity = 32; + int gridDensity = 50; - //Set pass band to stop at about 80% of channel bandwidth as starting point -// int passBandStop = (int)(channelBandwidth * 0.95); - int passBandStop = channelBandwidth; + int rolloffFrequency = (int)(alpha * (double)channelBandwidth); + + int passBandStop = channelBandwidth - rolloffFrequency; //Set the stop band to start at 40% into the adjacent channel, since we're oversampling by 2. - int stopBandStart = (int)(channelBandwidth * 1.40); + int stopBandStart = channelBandwidth + rolloffFrequency; - double passRipple = 0.1d; - double stopRipple = 0.01d; + double passRipple = 0.01d; + double stopRipple = 0.001d; //Approximately 90 dB of attenuation -// int order = FIRFilterSpecification.estimateBandPassOrder(sampleRate, passBandStart, passBandStop, -// passRipple, stopRipple); - int order = FIRFilterSpecification.estimateFilterOrder( sampleRate, passBandStop, - stopBandStart, 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) @@ -56,66 +54,39 @@ public class PolyphaseChannelizerFilterFactory order++; } + mLog.info("Filter Order: " + order); + FIRFilterSpecification specification = new PolyphaseChannelizerDesigner(order, gridDensity); - FIRFilterSpecification.FrequencyBand passBand = new FIRFilterSpecification.FrequencyBand(sampleRate, 0, - passBandStop, 1.0, passRipple); + FIRFilterSpecification.FrequencyBand passBand = new FIRFilterSpecification.FrequencyBand(sampleRate, passBandStart, + passBandStop, 1.0, passRipple, 1.0); + + //Use the filter order as the weighting for the pass band edge frequency + double weight = order; + + FIRFilterSpecification.FrequencyBand transitionBand = new FIRFilterSpecification.FrequencyBand(sampleRate, + channelBandwidth, channelBandwidth, OBJECTIVE_BAND_EDGE_COEFFICIENT_AMPLITUDE, passRipple, weight); + FIRFilterSpecification.FrequencyBand stopBand = new FIRFilterSpecification.FrequencyBand(sampleRate, stopBandStart, (int)(sampleRate / 2), 0.0, stopRipple); specification.addFrequencyBand(passBand); + specification.addFrequencyBand(transitionBand); 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) + try { - mLog.debug("Step Size: " + stepSize); - - RemezFIRFilterDesigner2 designer = new RemezFIRFilterDesigner2(specification); - - complete = true; //temporary - + RemezFIRFilterDesignerWithLagrange designer = new RemezFIRFilterDesignerWithLagrange(specification); + filter = designer.getImpulseResponse(); 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); - } + } + catch(Exception e) + { + mLog.error("Error designing filter", e); } return filter; @@ -125,7 +96,7 @@ public class PolyphaseChannelizerFilterFactory { mLog.debug("Starting ..."); - float[] filter = getFilter(2500000, 12500); + float[] filter = getFilter(200000, 12500, 0.21); mLog.debug("Finished ..."); } diff --git a/src/dsp/filter/fir/remez/RemezFIRFilterDesigner.java b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner.java index e5ca196b..dba1fd76 100644 --- a/src/dsp/filter/fir/remez/RemezFIRFilterDesigner.java +++ b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner.java @@ -212,15 +212,12 @@ public class RemezFIRFilterDesigner for( int k = 0; k < mExtremalIndices.size() - 1; k++ ) { - double cosineDelta = cosineOfFrequency - - mGrid.getCosineFrequencyGrid()[ mExtremalIndices.get( 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 this frequency is close to one of the polynomial points, use the polynomial point for the response if( Math.abs( cosineDelta ) < 1.0e-7 ) { - numerator = mIdealFrequencyResponse[ k ]; - denominator = 1.0; - break; + return mIdealFrequencyResponse[ k ]; } else { diff --git a/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java index 2890d939..521b4f65 100644 --- a/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java +++ b/src/dsp/filter/fir/remez/RemezFIRFilterDesigner2.java @@ -36,6 +36,7 @@ public class RemezFIRFilterDesigner2 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; + public static final double FREQUENCY_PROXIMITY_THRESHOLD = 1.0e-7; private FIRFilterSpecification mSpecification; private Grid mGrid; @@ -234,25 +235,27 @@ public class RemezFIRFilterDesigner2 for(int k = 0; k < mExtremalIndices.size() - 1; k++) { - BigDecimal cosineDelta = BigDecimal.valueOf(cosineOfFrequency - - mGrid.getCosineFrequencyGrid()[mExtremalIndices.get(k)]); + double cosineDelta = 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) + //If this frequency is close to one of the polynomial points, use the point as the response value + if(Math.abs(cosineDelta) < FREQUENCY_PROXIMITY_THRESHOLD) { - numerator = BigDecimal.valueOf(mIdealFrequencyResponse[k]); - denominator = BigDecimal.ONE; - break; + return mIdealFrequencyResponse[k]; } else { - BigDecimal dkOverCosineDelta = mD[k].divide(cosineDelta, 25, RoundingMode.HALF_UP); + BigDecimal dkOverCosineDelta = mD[k].divide(BigDecimal.valueOf(cosineDelta), 25, RoundingMode.HALF_UP); numerator = numerator.add(dkOverCosineDelta.multiply(BigDecimal.valueOf(mIdealFrequencyResponse[k]))); denominator = denominator.add(dkOverCosineDelta); } } + if(denominator.doubleValue() == 0.0) + { + int a = 0; + } + return numerator.divide(denominator, 25, RoundingMode.HALF_UP).doubleValue(); } diff --git a/src/dsp/filter/fir/remez/RemezFIRFilterDesignerWithLagrange.java b/src/dsp/filter/fir/remez/RemezFIRFilterDesignerWithLagrange.java new file mode 100644 index 00000000..d3fee4cf --- /dev/null +++ b/src/dsp/filter/fir/remez/RemezFIRFilterDesignerWithLagrange.java @@ -0,0 +1,686 @@ +/******************************************************************************* + * 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.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.math.BigDecimal; +import java.text.DecimalFormat; +import java.util.ArrayList; +import java.util.Iterator; +import java.util.List; + +public class RemezFIRFilterDesignerWithLagrange +{ + private final static Logger mLog = LoggerFactory.getLogger(RemezFIRFilterDesignerWithLagrange.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 mExtremalIndices; + + private double[] mD; + private double[] mGridErrors; + private double[] mGridFrequencyResponse; + private double[] mIdealFrequencyResponse; + private PolynomialFunctionLagrangeForm mLagrangeInterpolator; + + 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 RemezFIRFilterDesignerWithLagrange(FIRFilterSpecification specification) throws FilterDesignException + { + mSpecification = specification; + + design(); + } + + /** + * Designs the filter according to the specification + */ + private void design() throws FilterDesignException + { + mGrid = new Grid(mSpecification); + mExtremalIndices = getInitialExtremalIndices(); + + int iterationCount = 0; + + do + { + calculateGridFrequencyResponse(); + calculateGridError(); + findExtremalIndices(); + checkConvergence(); + + iterationCount++; + } + while(!mConverged && iterationCount < MAXIMUM_ITERATION_COUNT); + + if(mConverged) + { + //Update frequency response using the final set of extremal indices + calculateGridFrequencyResponse(); + } + else + { + throw new FilterDesignException("Unable to design filter after [" + iterationCount + "] iterations"); + } + } + + /** + * 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) + { + if(mLagrangeInterpolator == null) + { + throw new IllegalStateException("Lagrange interpolator is not setup - can't calculate frequency response"); + } + + return mLagrangeInterpolator.value(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 close to one of the polynomial points, use the polynomial point for the response +// if( Math.abs( cosineDelta ) < 1.0e-7 ) +// { +// return mIdealFrequencyResponse[ k ]; +// } +// 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 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() + { + BigDecimal[] b = calculateB(); + + calculateDelta(b); + + calculateC(); + + //Setup the lagrange polynomial so that we can determine the frequency response at the grid frequencies + //where mExtremalIndices are the polynomial points, the subset of cosine frequency grid identified by the + //extremal indices is the x value set and mIdealFrequencyResponse is the y axis value set + double[] xAxis = new double[mIdealFrequencyResponse.length]; + + for(int x = 0; x < mExtremalIndices.size() - 1; x++) + { + int index = mExtremalIndices.get(x); + + xAxis[x] = mGrid.getCosineFrequencyGrid()[index]; + } + + mLagrangeInterpolator = new PolynomialFunctionLagrangeForm(xAxis, mIdealFrequencyResponse); + + 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 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)); +// 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(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"); + } + } + + //Set the scale to the maximum scale we could see for a primitive decimal + 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 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 before checking for alternation theory - found:" + mExtremalIndices.size() + + " required:" + mSpecification.getExtremaCount()); + } + + //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 - found:" + mExtremalIndices.size() + " required:" + mSpecification.getExtremaCount()); + } + } + + /** + * 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); + + mLog.debug("Convergence: " + convergence); + + 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; + } +}