diff --git a/src/dsp/filter/FilterFactory.java b/src/dsp/filter/FilterFactory.java index 0889cd68..f45299d5 100644 --- a/src/dsp/filter/FilterFactory.java +++ b/src/dsp/filter/FilterFactory.java @@ -23,6 +23,7 @@ import dsp.filter.cic.ComplexPrimeCICDecimate; import dsp.filter.design.FilterDesignException; import dsp.filter.fir.FIRFilterSpecification; import dsp.filter.fir.remez.RemezFIRFilterDesigner; +import dsp.filter.fir.remez.RemezFIRFilterDesignerWithLagrange; import org.jtransforms.fft.FloatFFT_1D; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -36,7 +37,7 @@ import java.util.TreeSet; public class FilterFactory { private final static Logger mLog = LoggerFactory.getLogger(FilterFactory.class); - private static final double PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE = -6.020599842071533; //decibel(0.5) + private static final double PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE = -6.020599842071533; //decibel(0.5, 0.0) /** * Generates coefficients for a unity-gain, windowed low-pass filter @@ -695,6 +696,34 @@ public class FilterFactory return (float)(10.0 * Math.log10(Math.pow(real, 2.0) + Math.pow(imag, 2.0))); } + + public static float[] getRemezChannelizer(int channelBandwidth, int channels, int tapsPerChannel, double alpha, + double passRipple, double stopRipple) throws FilterDesignException + { + FIRFilterSpecification specification = FIRFilterSpecification.channelizerBuilder() + .sampleRate(channels * channelBandwidth * 4) + .channels(channels) + .channelBandwidth(channelBandwidth) + .tapsPerChannel(tapsPerChannel) + .alpha(alpha) + .passRipple(passRipple) + .stopRipple(stopRipple) + .build(); + + RemezFIRFilterDesignerWithLagrange designer = new RemezFIRFilterDesignerWithLagrange(specification); + + float[] taps = designer.getImpulseResponse(); + + double bandEdgeFrequency = (double)channelBandwidth / (double)(channels * channelBandwidth * 2); + + double response = decibel(designer.getFrequencyResponse(Math.cos(bandEdgeFrequency * Math.PI)), 0.0); + mLog.debug("Frequency Response at 1.0: " + response); + response = decibel(designer.getFrequencyResponse(Math.cos(bandEdgeFrequency * Math.PI * 2.0)), 0.0); + mLog.debug("Frequency Response at 2.0: " + response); + + return taps; + } + /** * Creates a windowed-sync (Nyquist) M/2 prototype FIR filter for use with a Polyphase Channelizer/Synthesizer. * The filter is designed for x2 oversampling of each channel and the filter cutoff frequency is incrementally @@ -720,8 +749,8 @@ public class FilterFactory * * @throws FilterDesignException if a filter cannot be designed with a band edge attenuation of -6.02 dB */ - public static float[] getChannelizer(int channelBandwidth, int channels, int tapsPerChannel, WindowType windowType, - boolean logResults) throws FilterDesignException + public static float[] getSincChannelizer(int channelBandwidth, int channels, int tapsPerChannel, WindowType windowType, + boolean logResults) throws FilterDesignException { int currentTapsPerChannel = tapsPerChannel; int sampleRate = channelBandwidth * channels * 2; diff --git a/src/dsp/filter/channelizer/ChannelizerViewer.java b/src/dsp/filter/channelizer/ChannelizerViewer.java index 0519486a..6d726a54 100644 --- a/src/dsp/filter/channelizer/ChannelizerViewer.java +++ b/src/dsp/filter/channelizer/ChannelizerViewer.java @@ -202,7 +202,7 @@ public class ChannelizerViewer extends JFrame private float[] getFilter() throws FilterDesignException { - return FilterFactory.getChannelizer(CHANNEL_BANDWIDTH, mChannelCount, 19, + return FilterFactory.getSincChannelizer(CHANNEL_BANDWIDTH, mChannelCount, 19, Window.WindowType.BLACKMAN_HARRIS_7, true); } diff --git a/src/dsp/filter/design/FilterViewer.java b/src/dsp/filter/design/FilterViewer.java index a7b973ee..14db2b9c 100644 --- a/src/dsp/filter/design/FilterViewer.java +++ b/src/dsp/filter/design/FilterViewer.java @@ -2,20 +2,12 @@ package dsp.filter.design; import dsp.filter.FilterFactory; import dsp.filter.Window; -import dsp.filter.fir.FIRFilterSpecification; -import dsp.filter.fir.remez.PolyphaseChannelizerFilterFactory; -import dsp.filter.fir.remez.RemezFIRFilterDesigner; -import dsp.filter.fir.remez.RemezFIRFilterDesigner2; -import dsp.filter.fir.remez.RemezFIRFilterDesignerWithLagrange; -import filter.Filter; import javafx.application.Application; import javafx.scene.Scene; import javafx.stage.Stage; import org.slf4j.Logger; import org.slf4j.LoggerFactory; -import java.util.Arrays; - public class FilterViewer extends Application { private final static Logger mLog = LoggerFactory.getLogger(FilterViewer.class); @@ -102,8 +94,10 @@ public class FilterViewer extends Application // taps = FilterFactory.getSinc(0.1, 51, Window.WindowType.BLACKMAN_HARRIS_7); try { - taps = FilterFactory.getChannelizer(12500, 50, 14, - Window.WindowType.BLACKMAN_HARRIS_7, true); +// taps = FilterFactory.getSincChannelizer(12500, 50, 14, +// Window.WindowType.BLACKMAN_HARRIS_7, true); + taps = FilterFactory.getRemezChannelizer(12500, 4, 22, 0.4, + 0.01, 0.0001); } catch(Exception e) { diff --git a/src/dsp/filter/fir/FIRFilterSpecification.java b/src/dsp/filter/fir/FIRFilterSpecification.java index da4aa40a..2ef0d820 100644 --- a/src/dsp/filter/fir/FIRFilterSpecification.java +++ b/src/dsp/filter/fir/FIRFilterSpecification.java @@ -18,7 +18,11 @@ ******************************************************************************/ package dsp.filter.fir; +import dsp.filter.FilterFactory; import dsp.filter.fir.remez.FIRLinearPhaseFilterType; +import org.apache.commons.lang3.Validate; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; import java.text.DecimalFormat; import java.util.ArrayList; @@ -27,6 +31,9 @@ import java.util.List; public class FIRFilterSpecification { + private final static Logger mLog = LoggerFactory.getLogger(FIRFilterSpecification.class); + + private static final double PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE = -6.020599842071533; //0.5 magnitude private static final DecimalFormat DECIMAL_FORMATTER = new DecimalFormat("0.00000"); private FIRLinearPhaseFilterType mRemezFilterType; @@ -143,6 +150,15 @@ public class FIRFilterSpecification return new BandPassBuilder(); } + /** + * Creates a Polyphase Channelizer filter specification builder that allows you to define the filter parameters and + * create a filter specification + */ + public static ChannelizerBuilder channelizerBuilder() + { + return new ChannelizerBuilder(); + } + public FIRLinearPhaseFilterType getFilterType() { return mRemezFilterType; @@ -693,57 +709,64 @@ public class FIRFilterSpecification /** * Builder for creating polyphase channelizer low-pass prototype filter. */ - public static class PolyphaseChannelizerBuilder + public static class ChannelizerBuilder { - private int mOrder; private int mGridDensity = 16; private int mSampleRate; + private int mChannelCount; private int mChannelBandwidth; - private int mPassFrequencyBegin = 0; - private int mPassFrequencyEnd; - private int mStopFrequencyBegin; + private int mTapsPerChannel; private double mStopRipple = 0.001; private double mPassRipple = 0.01; - private double mStopAmplitude = 0.0; - private double mPassAmplitude = 1.0; + private double mAlpha = 0.2; - public PolyphaseChannelizerBuilder() + public ChannelizerBuilder() { } /** * Sets the filter sample rate in hertz */ - public PolyphaseChannelizerBuilder sampleRate(int sampleRate) + public ChannelizerBuilder sampleRate(int sampleRate) { mSampleRate = sampleRate; return this; } + /** + * Number of Channels + * @param channels - must be a multiple of 2 + * @return + */ + public ChannelizerBuilder channels(int channels) + { + Validate.isTrue(channels % 2 == 0, "Channel count must be a multiple of 2"); + mChannelCount = channels; + return this; + } + /** * Sets the channel bandwidth in hertz */ - public PolyphaseChannelizerBuilder channelBandwidth(int channelBandwidth) + public ChannelizerBuilder 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. + * Number of filter taps per channel. Note: total filter length will be channels * tapsPerChannel */ - public PolyphaseChannelizerBuilder order(int order) + public ChannelizerBuilder tapsPerChannel(int tapsPerChannel) { - mOrder = order; + mTapsPerChannel = tapsPerChannel; return this; } /** * Sets the frequency grid density. Default is 16. */ - public PolyphaseChannelizerBuilder gridDensity(int density) + public ChannelizerBuilder gridDensity(int density) { mGridDensity = density; return this; @@ -752,7 +775,7 @@ public class FIRFilterSpecification /** * Specifies the ripple in dB for both stop bands */ - public PolyphaseChannelizerBuilder stopRipple(double rippleDb) + public ChannelizerBuilder stopRipple(double rippleDb) { mStopRipple = rippleDb; return this; @@ -761,66 +784,69 @@ public class FIRFilterSpecification /** * Specifies the ripple in dB for the pass band */ - public PolyphaseChannelizerBuilder passRipple(double rippleDb) + public ChannelizerBuilder passRipple(double rippleDb) { mPassRipple = rippleDb; return this; } /** - * Specifies the amplitude of the stop bands. Default is 0.0 + * Rolloff factor for the start band and stop band. + * @param alpha - value in range 0.0 - 1.0 */ - public PolyphaseChannelizerBuilder stopAmplitude(double amplitude) + public ChannelizerBuilder alpha(double alpha) { - mStopAmplitude = amplitude; - return this; - } + Validate.isTrue(0.0 <= alpha && alpha <= 1.0); + mAlpha = alpha; - /** - * 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; + FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL; - if(mOrder < 10) + int passFrequencyEnd = (int)(mChannelBandwidth * (1.0 - mAlpha)); + int stopFrequencyBegin = (int)(mChannelBandwidth * (1.0 + mAlpha)); + + int estimatedOrder = estimateFilterOrder(mSampleRate, passFrequencyEnd, stopFrequencyBegin, mPassRipple, + mStopRipple); + + // Ensure odd order since we're designing a Type 2 filter + if(estimatedOrder % 2 == 0) { - mOrder = estimateBandPassOrder(mSampleRate, mPassFrequencyBegin, mPassFrequencyEnd, - mPassRipple, mStopRipple); + estimatedOrder++; } - // Ensure even order since we're designing a Type 1 filter - if(mOrder % 2 == 1) - { - mOrder++; - } + int order = mChannelCount * mTapsPerChannel - 1; - FIRFilterSpecification spec = new FIRFilterSpecification(type, mOrder, mGridDensity); + mLog.debug("Order:" + order + " Sample Rate:" + mSampleRate + " Estimated Order:" + estimatedOrder + " Requested Order:" + order); - //Set pass band to stop at about 80% of channel bandwidth as starting point - int passBandStop = (int)(mChannelBandwidth * .80); + FIRFilterSpecification spec = new FIRFilterSpecification(type, order, mGridDensity); - //Set the stop band to start at 40% into the adjacent channel, since we're oversampling by 2. - int stopBandStart = (int)(mChannelBandwidth * 1.40); + double bandEdge = (double)mChannelBandwidth / (double)mSampleRate; + double passEnd = bandEdge * (1.0 - mAlpha); + double stopStart = bandEdge * (1.0 + mAlpha); - FrequencyBand passBand = new FrequencyBand(mSampleRate, 0, passBandStop, 1.0, mPassRipple); - FrequencyBand stopBand = new FrequencyBand(mSampleRate, stopBandStart, (int)(mSampleRate / 2), - 0.0, mStopRipple); + mLog.debug("Pass End:" + passEnd + " Band Edge:" + bandEdge + " Stop Begin:" + stopStart); + mLog.debug("Pass End:" + (int)(mSampleRate * passEnd) + + " Band Edge:" + (int)(mSampleRate * bandEdge) + + " Stop Begin:" + (int)(mSampleRate * stopStart)); + + FrequencyBand passBand = new FrequencyBand(0.0, passEnd, 1.0, mPassRipple); + FrequencyBand edgeBand = new FrequencyBand(bandEdge, bandEdge, 0.5, mPassRipple); + edgeBand.setWeight(8.0 * edgeBand.getWeight(mPassRipple)); + FrequencyBand stopBand = new FrequencyBand(stopStart, 0.5, 0.0, mStopRipple); spec.addFrequencyBand(passBand); + spec.addFrequencyBand(edgeBand); spec.addFrequencyBand(stopBand); + mLog.debug(spec.toString()); return spec; } } - // /** // * Converts the set of band amplitude values into a set of weights // * normalized to the largest band amplitude value @@ -992,6 +1018,14 @@ public class FIRFilterSpecification mWeight = weight; } + /** + * Explicitly sets the weight of this frequency band. + */ + public void setWeight(double weight) + { + mWeight = weight; + } + public int getGridSize() { return mGridSize; @@ -1018,6 +1052,8 @@ public class FIRFilterSpecification sb.append(DECIMAL_FORMATTER.format(getRippleAmplitude())); sb.append(" Grid Size: "); sb.append(mGridSize); + sb.append(" Weight: "); + sb.append(getWeight(mRippleDB)); return sb.toString(); }