Merges from master branch.

This commit is contained in:
Denny
2017-12-27 10:10:51 -05:00
parent 667c50174d
commit 236df8f3f4
262 changed files with 7701 additions and 49122 deletions

View File

@@ -26,18 +26,21 @@ import io.github.dsheirer.dsp.filter.cic.ComplexPrimeCICDecimate;
import io.github.dsheirer.dsp.filter.design.FilterDesignException;
import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification;
import io.github.dsheirer.dsp.filter.fir.remez.RemezFIRFilterDesigner;
import io.github.dsheirer.dsp.filter.fir.remez.RemezFIRFilterDesignerWithLagrange;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.Arrays;
import java.util.Set;
import java.util.TreeSet;
/**
* Utility methods for designing various types of filters.
*/
public class FilterFactory
{
private final static Logger mLog =
LoggerFactory.getLogger(FilterFactory.class);
private final static Logger mLog = LoggerFactory.getLogger(FilterFactory.class);
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
@@ -48,43 +51,49 @@ public class FilterFactory
* @param window - to apply against the coefficients
* @return
*/
public static float[] getSinc(int sampleRate,
long frequency,
int length,
Window.WindowType window)
public static float[] getSinc(int sampleRate, long frequency, int length, Window.WindowType window)
{
//Ensure we have an odd length
assert (length % 2 == 0);
int evenLength = length % 2 == 0 ? length : length + 1;
//Get unity response array (one element longer to align with IDFT size)
float[] frequencyResponse = getUnityResponseArray(sampleRate, frequency, length + 1);
float[] frequencyResponse = getUnityResponseArray(sampleRate, frequency, evenLength);
//Apply Inverse DFT against frequency response unity values, leaving the
//IDFT bin results in the frequency response array
FloatFFT_1D idft = new FloatFFT_1D(length + 1);
idft.realInverseFull(frequencyResponse, false);
FloatFFT_1D idft = new FloatFFT_1D(evenLength);
idft.realInverseFull(frequencyResponse, true);
//Transfer the IDFT results to the return array
float[] coefficients = new float[length];
int middleCoefficient = (int) (length / 2);
int middleCoefficient = (int)(length / 2);
//Bin 0 of the idft is our center coefficient
coefficients[middleCoefficient] = frequencyResponse[0];
//The remaining idft bins from 1 to (middle - 1) are the mirror image
//coefficients
for(int x = 1; x < middleCoefficient; x++)
if(length % 2 == 0) //Even length
{
coefficients[middleCoefficient + x] = frequencyResponse[2 * x];
coefficients[middleCoefficient - x] = frequencyResponse[2 * x];
//The remaining idft bins from 1 to (middle - 1) are the mirror image
//coefficients
for(int x = 0; x < middleCoefficient; x++)
{
coefficients[middleCoefficient + x] = frequencyResponse[2 * x];
coefficients[middleCoefficient - x] = frequencyResponse[2 * x];
}
}
else
{
//Bin 0 of the idft is our center coefficient
coefficients[middleCoefficient] = frequencyResponse[0];
//The remaining idft bins from 1 to (middle - 1) are the mirror image
//coefficients
for(int x = 1; x < middleCoefficient; x++)
{
coefficients[middleCoefficient + x] = frequencyResponse[2 * x];
coefficients[middleCoefficient - x] = frequencyResponse[2 * x];
}
}
//Apply the window against the coefficients
coefficients = Window.apply(window, coefficients);
//Normalize to unity (1) gain
coefficients = normalize(coefficients);
return coefficients;
}
@@ -95,7 +104,7 @@ public class FilterFactory
* @param coefficients
* @return
*/
private static float[] normalize(float[] coefficients)
public static float[] normalize(float[] coefficients)
{
float accumulator = 0;
@@ -128,17 +137,15 @@ public class FilterFactory
* @param length
* @return
*/
public static float[] getUnityResponseArray(int sampleRate,
long frequency,
int length)
public static float[] getUnityResponseArray(int sampleRate, long frequency, int length)
{
float[] unityArray = new float[length * 2];
int binCount = (int) ((Math.round(
(float) frequency / (float) sampleRate * (float) length)));
if(length % 2 == 0) //even length
{
int binCount = (int)((Math.round((float)frequency / (float)sampleRate * (float)length)));
for(int x = 0; x < binCount; x++)
{
unityArray[x] = 1.0f;
@@ -147,9 +154,11 @@ public class FilterFactory
}
else //odd length
{
int binCount = (int)((Math.round((float)frequency / (float)sampleRate * (float)length + 0.5)));
unityArray[0] = 1.0f;
for(int x = 1; x <= binCount; x++)
for(int x = 1; x < binCount; x++)
{
unityArray[x] = 1.0f;
unityArray[length - x] = 1.0f;
@@ -159,25 +168,6 @@ public class FilterFactory
return unityArray;
}
public static float[] getSine(float sampleRate, float frequency, int length)
{
float[] retVal = new float[length];
float radianFrequency = (float) (2.0d * Math.PI * (frequency / sampleRate));
int middle = (int) (length / 2);
for(int x = 0; x < middle; x++)
{
float val = (float) Math.sin(radianFrequency * x);
retVal[middle + x] = val;
retVal[middle - x] = -val;
}
return retVal;
}
/**
* Applies a repeating sequence of 1, -1 to the coefficients to invert
* the frequency response of the filter. Used in converting a low-pass
@@ -211,25 +201,9 @@ public class FilterFactory
* @param windowType - window to apply against the generated coefficients
* @return
*/
public static float[] getLowPass(int sampleRate,
long cutoff,
int filterLength,
Window.WindowType windowType)
public static float[] getLowPass(int sampleRate, long cutoff, int filterLength, Window.WindowType windowType)
{
if(filterLength % 2 == 0) //even length
{
float[] values = getSinc(sampleRate, cutoff, filterLength + 2, windowType);
//throw away the 0 index and the last index
return Arrays.copyOfRange(values, 1, values.length - 2);
}
else
{
float[] values = getSinc(sampleRate, cutoff, filterLength + 1, windowType);
//throw away the 0 index
return Arrays.copyOfRange(values, 1, values.length);
}
return getSinc(sampleRate, cutoff, filterLength, windowType);
}
/**
@@ -246,26 +220,17 @@ public class FilterFactory
* - passFrequency < stopFrequency
* - stopFrequency <= sampleRate/2
*/
public static float[] getLowPass(int sampleRate,
int passFrequency,
int stopFrequency,
int attenuation,
Window.WindowType windowType,
boolean forceOddLength)
public static float[] getLowPass(int sampleRate, int passFrequency, int stopFrequency, int attenuation,
Window.WindowType windowType, boolean forceOddLength)
{
if(stopFrequency < passFrequency || stopFrequency > (sampleRate / 2))
{
throw new IllegalArgumentException("FilterFactory - low pass "
+ "filter pass frequency [" + passFrequency + "] must be "
+ "less than the stop frequency [" + stopFrequency + "] "
+ "and must be less than half [" + (int) (sampleRate / 2) +
"] of the sample rate [" + sampleRate + "]");
throw new IllegalArgumentException("FilterFactory - low pass filter pass frequency [" + passFrequency +
"] must be less than the stop frequency [" + stopFrequency + "] and must be less than half [" +
(int)(sampleRate / 2) + "] of the sample rate [" + sampleRate + "]");
}
int tapCount = getTapCount(sampleRate,
passFrequency,
stopFrequency,
attenuation);
int tapCount = getTapCount(sampleRate, passFrequency, stopFrequency, attenuation);
if(forceOddLength)
{
@@ -289,34 +254,21 @@ public class FilterFactory
* @param windowType - window to apply against the generated coefficients
* @return
*/
public static float[] getHighPass(int sampleRate,
long cutoff,
int filterLength,
Window.WindowType windowType)
public static float[] getHighPass(int sampleRate, long cutoff, int filterLength, Window.WindowType windowType)
{
//Convert the high frequency cutoff to its low frequency cutoff
//equivalent, so that when we generate the low pass filter, prior to
//inversion, its at the correct frequency
long convertedCutoff = sampleRate / 2 - cutoff;
return invert(getSinc(sampleRate,
convertedCutoff,
filterLength,
windowType));
return invert(getSinc(sampleRate, convertedCutoff, filterLength, windowType));
}
public static float[] getHighPass(int sampleRate,
long stopFrequency,
long passFrequency,
int attenuation,
Window.WindowType windowType,
boolean forceOddLength)
public static float[] getHighPass(int sampleRate, long stopFrequency, long passFrequency, int attenuation,
Window.WindowType windowType, boolean forceOddLength)
{
/* reverse the stop and pass frequency to get the low pass variant */
int tapCount = getTapCount(sampleRate,
stopFrequency,
passFrequency,
attenuation);
int tapCount = getTapCount(sampleRate, stopFrequency, passFrequency, attenuation);
if(forceOddLength)
{
@@ -360,11 +312,10 @@ public class FilterFactory
}
/**
* Determines the number of fir filter taps required to produce the
* specified frequency response with passband ripple near .1dB.
* Determines the number of fir filter taps required to produce the specified frequency response with passband
* ripple near 0.1dB.
*
* Implements the algorithm from Understanding Digital Signal Processing, 3e
* , Lyons, section 5.10.5.
* Implements the algorithm from Understanding Digital Signal Processing, 3e, Lyons, section 5.10.5.
*
* @param sampleRate in hertz
* @param pass pass frequency in hertz
@@ -372,21 +323,17 @@ public class FilterFactory
* @param attenuation in dB
* @return
*/
public static int getTapCount(int sampleRate,
long pass,
long stop,
int attenuation)
public static int getTapCount(int sampleRate, long pass, long stop, int attenuation)
{
double frequency = ((double) stop - (double) pass) / (double) sampleRate;
double frequency = ((double)stop - (double)pass) / (double)sampleRate;
return (int) (Math.round((double) attenuation / (22.0d * frequency)));
return (int)(Math.round((double)attenuation / (22.0d * frequency)));
}
public static ComplexPrimeCICDecimate getDecimationFilter(int sampleRate,
int decimatedRate, int order, int passFrequency, int attenuation,
Window.WindowType windowType)
public static ComplexPrimeCICDecimate getDecimationFilter(int sampleRate, int decimatedRate, int order,
int passFrequency, int attenuation, Window.WindowType windowType)
{
int decimationRate = (int) (sampleRate / decimatedRate);
int decimationRate = (int)(sampleRate / decimatedRate);
return new ComplexPrimeCICDecimate(decimationRate, order,
passFrequency, attenuation, windowType);
@@ -401,10 +348,7 @@ public class FilterFactory
* polyphase filter decimation chain
* @throws - AssertionException if sample rate is not a multiple of 48 kHz
*/
public static int[] getPolyphaseDecimationRates(int sampleRate,
int decimatedRate,
long passFrequency,
long stopFrequency)
public static int[] getPolyphaseDecimationRates(int sampleRate, int decimatedRate, long passFrequency, long stopFrequency)
{
int[] rates;
@@ -414,7 +358,7 @@ public class FilterFactory
+ "integer multiple of sample rate");
}
int decimation = (int) (sampleRate / decimatedRate);
int decimation = (int)(sampleRate / decimatedRate);
//Decimation rates below 20 will use single stage polyphase filter
if(decimation < 20)
@@ -432,11 +376,6 @@ public class FilterFactory
int stage1 = findClosest(optimalStage1, factors);
// mLog.info( "Decimation rate [" + decimation +
// "] stage1 optimal [" + optimalStage1 +
// "] stage1 actual [" + stage1 +
// "]");
if(stage1 == decimation || stage1 == 1)
{
rates = new int[1];
@@ -447,7 +386,7 @@ public class FilterFactory
{
rates = new int[2];
int stage2 = (int) (decimation / stage1);
int stage2 = (int)(decimation / stage1);
if(stage1 > stage2)
{
@@ -465,8 +404,7 @@ public class FilterFactory
}
/**
* Finds the factor that is closest to the desired factor, from an ordered
* list of factors.
* Finds the factor that is closest to the desired factor, from an ordered list of factors.
*/
private static int findClosest(int desiredFactor, Set<Integer> factors)
{
@@ -488,9 +426,8 @@ public class FilterFactory
}
/**
* Determines the factors that make up an integer. Uses a brute force
* method to iterate all integers from 1 to value, determining which factors
* are evenly divisible into the value.
* Determines the factors that make up an integer. Uses a brute force method to iterate all integers from 1 to
* value, determining which factors are evenly divisible into the value.
*
* @param value - integer decimation value
* @return - ordered set of factors for value
@@ -498,11 +435,11 @@ public class FilterFactory
private static Set<Integer> getFactors(int value)
{
Set<Integer> factors = new TreeSet<Integer>();
/* Brute force */
for(int x = 1; x <= value; x++)
{
int remainder = (int) (value / x);
int remainder = (int)(value / x);
if(remainder * x == value)
{
@@ -514,14 +451,10 @@ public class FilterFactory
}
/**
* Determines the optimal decimation rate for stage 1 of a two-stage
* poly-phase decimation filter chain, to produce a final sample rate of
* 48 kHz using a pass bandwidth of 25 kHz.
* Determines the optimal decimation rate for stage 1 of a two-stage poly-phase decimation filter chain, to produce
* a final sample rate of 48 kHz using a pass bandwidth of 25 kHz. Use for total decimation rates of 20 or higher.
*
* Use for total decimation rates of 20 or higher.
*
* Implements the algorithm described in Lyons, Understanding Digital Signal
* Processing, 3e, section 10.2.1, page 511.
* Implements the algorithm described in Lyons, Understanding Digital Signal Processing, 3e, section 10.2.1, page 511.
*
* @param sampleRate
* @param decimation
@@ -529,42 +462,29 @@ public class FilterFactory
* @return optimum integer decimation rate for the first stage decimation
* filter
*/
public static int getOptimalStageOneRate(int sampleRate,
int decimation,
long passFrequency,
long stopFrequency)
public static int getOptimalStageOneRate(int sampleRate, int decimation, long passFrequency, long stopFrequency)
{
double ratio = getBandwidthRatio(passFrequency, stopFrequency);
double numerator = 1.0d - (
Math.sqrt((double) decimation * ratio / (2.0d - ratio)));
double numerator = 1.0d - (Math.sqrt((double)decimation * ratio / (2.0d - ratio)));
double denominator = 2.0d - (ratio * (decimation + 1.0d));
int retVal = (int) (2.0d * decimation * (numerator / denominator));
// mLog.info( "Optimal Stage 1 Decimation - rate [" + sampleRate +
// "] pass [" + passFrequency +
// "] bw ratio [" + ratio +
// "] optimal [" + retVal +
// "]" );
int retVal = (int)(2.0d * decimation * (numerator / denominator));
return retVal;
}
/**
* Determines the F ratio as described in Lyons, Understanding Digital
* Signal Processing, 3e, section 10.2.1, page 511
*
* Used in conjunction with the optimal stage one decimation rate method
* above.
* Determines the F ratio as described in Lyons, Understanding Digital Signal Processing, 3e, section 10.2.1, page 511
* Used in conjunction with the optimal stage one decimation rate method above.
*/
private static double getBandwidthRatio(long passFrequency, long stopFrequency)
{
assert (passFrequency < stopFrequency);
return ((double) (stopFrequency - passFrequency) /
(double) stopFrequency);
return ((double)(stopFrequency - passFrequency) /
(double)stopFrequency);
}
/**
@@ -575,14 +495,11 @@ public class FilterFactory
* @param outputSampleRate
* @return
*/
public static float[] getCICCleanupFilter(int outputSampleRate,
int passFrequency,
int attenuation,
Window.WindowType window)
public static float[] getCICCleanupFilter(int outputSampleRate, int passFrequency, int attenuation, Window.WindowType window)
{
int taps = getTapCount(outputSampleRate, passFrequency, passFrequency + 1500,
attenuation);
/* Make tap count odd */
if(taps % 2 == 0)
{
@@ -599,7 +516,7 @@ public class FilterFactory
//Transfer the IDFT results to the odd length return array
float[] coefficients = new float[taps];
int middleCoefficient = (int) (taps / 2);
int middleCoefficient = (int)(taps / 2);
//Bin 0 of the idft is our center coefficient
coefficients[middleCoefficient] = frequencyResponse[0];
@@ -620,26 +537,22 @@ public class FilterFactory
return coefficients;
}
public static float[] getCICResponseArray(int sampleRate,
int frequency,
int length)
public static float[] getCICResponseArray(int sampleRate, int frequency, int length)
{
float[] cicArray = new float[length * 2];
int binCount = (int) ((Math.round(
(double) frequency / (double) sampleRate * 2.0d * (double) length)));
int binCount = (int)((Math.round(
(double)frequency / (double)sampleRate * 2.0d * (double)length)));
cicArray[0] = 1.0f;
float unityResponse = (float) (Math.sin(1.0d / (double) length) /
(1.0d / (double) length));
float unityResponse = (float)(Math.sin(1.0d / (double)length) /
(1.0d / (double)length));
for(int x = 1; x <= binCount; x++)
{
/* Calculate unity response plus amplification for droop */
float compensated = 1.0f + (unityResponse -
(float) (Math.sin((double) x / (double) length) /
((double) x / (double) length)));
/* Calculate unity response plus amplification for droop */
float compensated = 1.0f + (unityResponse - (float)(Math.sin((double)x / (double)length) / ((double)x / (double)length)));
cicArray[x] = compensated;
cicArray[length - x] = compensated;
@@ -649,28 +562,29 @@ public class FilterFactory
}
/**
* Creates root raised cosine filter coefficients with a tap count equal
* to the symbols x samplesPerSymbol + 1.
* Creates root raised cosine filter coefficients with a tap count equal to the symbols x samplesPerSymbol + 1.
*
* Symbol count should be an even integer.
*
* Ported to java from gnuradio/filter/firdes.cc
*
* For 40db attenuation, calculate the number of symbols based on the
* following formula:
* For 40db attenuation, calculate the number of symbols based on the following formula:
* Symbols = -44 * alpha + 33
*
* Symbols = -44 * alpha + 33
* Polyphase Channelizer notes:
* -Set samples Per Symbol at 2 (or more) * number of channels
* -Set symbolCount to sufficient size to produce required attenuation
* -Set symbol rate in hertz
* -Set alpha = (desired channel bandwidth / (symbol rate * samples per symbol) - 1.0
* -Channel bandwidth must be larger than symbol rate * samples per symbol, in order for there to be a positive alpha
*
* @param samplesPerSymbol - number of samples per symbol
* @param symbols - number of symbols - must be even
* @param alpha - roll-off factor
* @param symbolCount - number of symbol intervals for the filter. This directly impacts the filter size
* @param alpha - roll-off factor.
* @return - filter coefficients
*/
public static float[] getRootRaisedCosine(int samplesPerSymbol,
int symbols,
float alpha)
public static float[] getRootRaisedCosine(int samplesPerSymbol, int symbolCount, float alpha)
{
int taps = samplesPerSymbol * symbols + 1;
int taps = samplesPerSymbol * symbolCount;
float scale = 0;
@@ -678,10 +592,10 @@ public class FilterFactory
for(int x = 0; x < taps; x++)
{
float index = (float) x - ((float) taps / 2.0f);
float index = (float)x - ((float)taps / 2.0f);
float x1 = (float) Math.PI * index / (float) samplesPerSymbol;
float x2 = 4.0f * alpha * index / (float) samplesPerSymbol;
float x1 = (float)Math.PI * index / (float)samplesPerSymbol;
float x2 = 4.0f * alpha * index / (float)samplesPerSymbol;
float x3 = x2 * x2 - 1.0f;
float numerator, denominator;
@@ -690,17 +604,15 @@ public class FilterFactory
{
if(x != taps / 2)
{
numerator = (float) Math.cos((1.0 + alpha) * x1) +
(float) Math.sin((1.0f - alpha) * x1) /
(4.0f * alpha * index / (float) samplesPerSymbol);
numerator = (float)Math.cos((1.0 + alpha) * x1) +
(float)Math.sin((1.0f - alpha) * x1) / (4.0f * alpha * index / (float)samplesPerSymbol);
}
else
{
numerator = (float) Math.cos((1.0f + alpha) * x1) +
(1.0f - alpha) * (float) Math.PI / (4.0f * alpha);
numerator = (float)Math.cos((1.0f + alpha) * x1) + (1.0f - alpha) * (float)Math.PI / (4.0f * alpha);
}
denominator = x3 * (float) Math.PI;
denominator = x3 * (float)Math.PI;
}
else
{
@@ -713,14 +625,11 @@ public class FilterFactory
x3 = (1.0f - alpha) * x1;
x2 = (1.0f + alpha) * x1;
numerator = (float) (Math.sin(x2) * (1.0f + alpha) * Math.PI -
Math.cos(x3) * ((1.0 - alpha) * Math.PI *
(double) samplesPerSymbol) / (4 * alpha * index) +
Math.sin(x3) * (double) samplesPerSymbol *
(double) samplesPerSymbol / (4.0 * alpha * index * index));
numerator = (float)(Math.sin(x2) * (1.0f + alpha) * Math.PI - Math.cos(x3) * ((1.0 - alpha) * Math.PI *
(double)samplesPerSymbol) / (4.0 * alpha * index) + Math.sin(x3) * (double)samplesPerSymbol *
(double)samplesPerSymbol / (4.0 * alpha * index * index));
denominator = (float) (-32.0 * Math.PI * alpha * alpha * index /
(double) samplesPerSymbol);
denominator = (float)(-32.0 * Math.PI * alpha * alpha * index / (double)samplesPerSymbol);
}
coefficients[x] = 4.0f * alpha * numerator / denominator;
@@ -740,10 +649,10 @@ public class FilterFactory
}
/**
* Creates a filter from the filter specification
* Creates a filter from the filter specification using the remez exchange design algorithm
*
* @param specification
* @return
* @return filter coefficients
* @throws FilterDesignException if the filter cannot be designed
*/
public static float[] getTaps(FIRFilterSpecification specification) throws FilterDesignException
@@ -757,4 +666,248 @@ public class FilterFactory
return null;
}
/**
* Evaluates the filter to determine the magnitude response (dB) at the specified frequency.
*
* @param filter to evaluate
* @param frequency to calculate (0 <> 0.5)
* @return magnitude frequency response in decibels
*/
public static double evaluate(float[] filter, double frequency)
{
double real = 0.0;
double imag = 0.0;
for(int x = 0; x < filter.length; x++)
{
real += filter[x] * Math.cos(Math.PI * frequency * (double)x);
imag += filter[x] * Math.sin(Math.PI * frequency * (double)x);
}
return decibel(real, imag);
}
/**
* Calculates the decibel magnitude of the real and imaginary complex sample
* @param real (x-axis) value
* @param imag (y-axis) value
* @return magnitude in decibels
*/
public static double decibel(double real, double imag)
{
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
* adjusted to achieve a -6.02 dB attenuation at the channel band edge to enable perfect reconstruction of adjacent
* channels.
*
* The odd length sync filter that meets the design objectives will be increased in length by pre-padding a zero
* coefficient to produce a filter with a length that is an integral of the channel count.
*
* Note: the tapsPerChannel value will be iteratively increased (up to 10x) if a filter cannot be designed to meet
* the channel band edge attenuation objective.
*
* Note: a channelizer filter designed using a Blackman-Harris-7 window and 18-21 taps per channel creates a filter
* with attenuation in range of -60 to -80 dB at the edge of the 2x sampled channel.
*
* @param channelBandwidth - full bandwidth of the channel in hertz
* @param channels - number of channels
* @param tapsPerChannel - desired filter tap count per channel (may be increased to meet objective)
* @param windowType - to use when designing the windowed-sync filter
* @param logResults - to log a debug summary of the filter characteristics
*
* @return filter coefficients
*
* @throws FilterDesignException if a filter cannot be designed with a band edge attenuation of -6.02 dB
*/
public static float[] getSincChannelizer(double channelBandwidth, int channels, int tapsPerChannel,
Window.WindowType windowType, boolean logResults) throws FilterDesignException
{
int currentTapsPerChannel = tapsPerChannel;
int sampleRate = (int)(channelBandwidth * (double)channels * 2.0);
int filterLength = (channels * currentTapsPerChannel) - 1;
double cutoffFrequency = channelBandwidth / (double)sampleRate;
double bandEdge = cutoffFrequency;
double increment = cutoffFrequency * 0.1;
//Get an initial filter and band edge frequency response using the channel bandwidth as a cutoff
float[] taps = FilterFactory.getSinc(cutoffFrequency, filterLength, windowType);
double response = FilterFactory.evaluate(taps, bandEdge);
//Set cutoff adjustment threshold - we'll test cutoff frequencies to around 1 hertz resolution
double incrementThreshold = 1.0 / (double)sampleRate;
//Iteratively evaluate filters to find the filter with the highest cutoff frequency that meets the objective
while(increment > incrementThreshold)
{
//If the current cutoff meets the objective, test if a higher cutoff also meets the objective
if(matchesObjective(response, PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE) &&
(cutoffFrequency + increment <= bandEdge))
{
float[] higherCutoffTaps = FilterFactory.getSinc(cutoffFrequency + increment, filterLength, windowType);
double higherCutoffResponse = FilterFactory.evaluate(higherCutoffTaps, bandEdge);
if(matchesObjective(higherCutoffResponse, PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE))
{
cutoffFrequency += increment;
taps = higherCutoffTaps;
response = higherCutoffResponse;
}
else
{
increment /= 2.0;
}
}
//If the current cutoff meets the objective, decrease the increment to test for smaller resolution cutoffs
else if(matchesObjective(response, PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE))
{
increment /= 2.0;
}
//Decrease the cutoff frequency - we can't meet the objective
else
{
cutoffFrequency -= increment;
if(cutoffFrequency <= 0)
{
if(logResults)
{
mLog.debug("Warning: cannot design channelizer filter with tap count [" + currentTapsPerChannel +
"] increasing tap count and starting over");
}
currentTapsPerChannel++;
//After we've attempted increasing the taps per channel 10x, give up
if(currentTapsPerChannel > (tapsPerChannel + 10))
{
throw new FilterDesignException("Couldn't design filter with taps per channel count in the " +
"range of " + tapsPerChannel + " - " + (tapsPerChannel + 10));
}
filterLength = channels * currentTapsPerChannel - 1;
cutoffFrequency = ((double)channelBandwidth / (double)sampleRate);
increment = cutoffFrequency * 0.1;
}
taps = FilterFactory.getSinc(cutoffFrequency, filterLength, windowType);
response = FilterFactory.evaluate(taps, bandEdge);
}
}
//This will probably never happen, but ...
if(!matchesObjective(response, PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE))
{
throw new FilterDesignException("Cannot design filter to specifications");
}
//This is an odd length filter - increase the length by 1 by pre-padding a zero coefficient
float[] extendedTaps = new float[taps.length + 1];
System.arraycopy(taps, 0, extendedTaps, 1, taps.length);
if(logResults)
{
mLog.debug("Polyphase Channelizer Filter Design Summary");
mLog.debug("-----------------------------------------------------");
mLog.debug("Input Sample Rate: " + sampleRate);
mLog.debug("Channel Bandwidth: " + channelBandwidth);
mLog.debug("Channels: " + channels);
mLog.debug("Window Type: " + windowType.name());
mLog.debug("Taps Per Channel - Requested:" + tapsPerChannel + " Actual:" + ((double)extendedTaps.length / (double)channels));
mLog.debug("Filter Length: " + (extendedTaps.length));
mLog.debug("Requested Cutoff Frequency: " + (sampleRate * bandEdge));
mLog.debug("Actual Cutoff Frequency: " + (sampleRate * cutoffFrequency));
mLog.debug("Attenuation at 1.0 OBJECTIVE: " + PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE);
mLog.debug("Attenuation at 1.00 Channels: " + evaluate(taps, bandEdge * 1.00));
mLog.debug("Attenuation at 1.25 Channels: " + evaluate(taps, bandEdge * 1.25));
mLog.debug("Attenuation at 1.50 Channels: " + evaluate(taps, bandEdge * 1.50));
mLog.debug("Attenuation at 1.75 Channels: " + evaluate(taps, bandEdge * 1.75));
mLog.debug("Attenuation at 2.00 Channels: " + evaluate(taps, bandEdge * 2.00));
}
return extendedTaps;
}
/**
* Creates a windowed-sync (Nyquist) filter.
*
* @param cutoff frequency (0.0 <> 0.5)
* @param length of the filter - must be odd-length
* @param windowType to use in designing the filter
* @return filter coefficients.
* @throws FilterDesignException if the requested length is not odd
*/
public static float[] getSinc(double cutoff, int length, Window.WindowType windowType) throws FilterDesignException
{
if(length % 2 == 0)
{
throw new FilterDesignException("Sinc filters must be odd-length");
}
float[] coefficients = new float[length];
int half = length / 2;
double[] window = Window.getWindow(windowType, length);
double scalor = 2.0 * cutoff;
double piScalor = Math.PI * scalor;
coefficients[half] = (float)(1.0 * scalor);
for(int x = 1; x <= half; x++)
{
double a = piScalor * x;
double coefficient = scalor * Math.sin(a) / a;
coefficient *= window[half + x];
coefficients[half + x] = (float)coefficient;
coefficients[half - x] = (float)coefficient;
}
return coefficients;
}
/**
* Compares two doubles for equals and avoids any rounding error that are present.
* @param a value to compare
* @param objective value to compare
* @return true if they are equivalent within the margin of error
*/
private static boolean matchesObjective(double a, double objective)
{
return Math.abs(a - objective) < 0.0000000001;
}
}

View File

@@ -1,234 +1,439 @@
/*******************************************************************************
* SDR Trunk
* Copyright (C) 2014 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.Arrays;
/**
* Window creator
*
* Blackman, Hamming and Hanning windows
* based on windowing methods described at:
* http://www.mstarlabs.com/dsp/goertzel/goertzel.html
* and
* Understanding Digital Signal Processing 3e, Lyons, p.91
* and
* Digital Filters for Everyone, Allred, p. 132
*
* Window factory and corresponding utility methods.
*/
public class Window
{
public static double[] getWindow( WindowType type, int length )
private final static Logger mLog = LoggerFactory.getLogger(Window.class);
private static final double TWO_PI = Math.PI * 2.0d;
private static final double FOUR_PI = Math.PI * 4.0d;
private static final double SIX_PI = Math.PI * 6.0d;
private static final double EIGHT_PI = Math.PI * 8.0d;
/**
* Creates a window of the specified type and length.
*
* @param type of window
* @param length of window
* @return window of the specified type and length
*/
public static double[] getWindow(WindowType type, int length)
{
switch( type )
{
case BLACKMAN:
return getBlackmanWindow( length );
case COSINE:
return getCosineWindow( length );
case HAMMING:
return getHammingWindow( length );
case HANNING:
return getHanningWindow( length );
case NONE:
default:
return getRectangularWindow( length );
}
switch(type)
{
case BLACKMAN:
return getBlackman(length);
case BLACKMAN_HARRIS_4:
return getBlackmanHarris4(length);
case BLACKMAN_HARRIS_7:
return getBlackmanHarris7(length);
case BLACKMAN_NUTALL:
return getBlackmanNutall(length);
case COSINE:
return getCosine(length);
case FLAT_TOP:
return getFlatTop(length);
case HAMMING:
return getHamming(length);
case HANN:
return getHann(length);
case NUTALL:
return getNutall(length);
case NONE:
default:
return getRectangular(length);
}
}
public static double[] getRectangularWindow( int length )
/**
* Creates an all-pass or rectangular window.
*
* @param length of the window
* @return rectangular window
*/
public static double[] getRectangular(int length)
{
double[] coefficients = new double[ length ];
double[] coefficients = new double[length];
Arrays.fill( coefficients, 1.0d );
return coefficients;
}
public static double[] getCosineWindow( int length )
{
double[] coefficients = new double[ length ];
Arrays.fill(coefficients, 1.0d);
if( length % 2 == 0 ) //Even length
{
int half = (int)( ( length - 1 ) / 2 );
for( int x = -half; x < length / 2 + 1; x++ )
{
coefficients[ x + half ] = Math.cos(
( (double)x * Math.PI ) / ( (double)length + 1.0d ) );
}
}
else //Odd length
{
int half = (int) length / 2;
for( int x = -half; x < half + 1; x++ )
{
coefficients[ x + half ] = Math.cos(
( (double)x * Math.PI ) / ( (double)length + 1.0d ) );
}
}
return coefficients;
return coefficients;
}
public static double[] getBlackmanWindow( int length )
/**
* Creates a cosine window of the specified length.
*
* @param length of the window
* @return cosine window
*/
public static double[] getCosine(int length)
{
double[] coefficients = new double[ length ];
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = .426591D -
( .496561D * Math.cos( ( Math.PI * 2.0D * (double)x ) / (double)( length - 1 ) ) +
( .076848D * Math.cos( ( Math.PI * 4.0D * (double)x ) / (double)( length - 1 ) ) ) );
}
double[] coefficients = new double[length];
return coefficients;
}
public static double[] getHammingWindow( int length )
{
double[] coefficients = new double[ length ];
if( length % 2 == 0 ) //Even length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = .54D - ( .46D *
Math.cos( ( Math.PI * ( 2.0D * (double)x + 1.0d ) ) / (double)( length - 1 ) ) );
}
}
else //Odd length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = .54D - ( .46D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) );
}
}
if(length % 2 == 0) //Even length
{
int half = (int)((length - 1) / 2);
return coefficients;
for(int x = -half; x < length / 2 + 1; x++)
{
coefficients[x + half] = Math.cos(
((double)x * Math.PI) / ((double)length + 1.0d));
}
}
else //Odd length
{
int half = (int)length / 2;
for(int x = -half; x < half + 1; x++)
{
coefficients[x + half] = Math.cos(
((double)x * Math.PI) / ((double)length + 1.0d));
}
}
return coefficients;
}
public static double[] getHanningWindow( int length )
/**
* Creates a 3-term Blackman window.
*
* @param length of the window
* @return Blackman window
*/
public static double[] getBlackman(int length)
{
double[] coefficients = new double[ length ];
double[] coefficients = new double[length];
if( length % 2 == 0 ) //Even length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = .5D - ( .5D *
Math.cos( ( Math.PI * ( 2.0D * (double)x + 1 ) ) / (double)( length - 1 ) ) );
}
}
else //Odd length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = .5D - ( .5D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) );
}
}
double denominator = length - 1;
double a0 = 0.426590713672;
double a1 = 0.496560619089;
double a2 = 0.0768486672399;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 -
(a1 * Math.cos((TWO_PI * (double)x) / denominator)) +
(a2 * Math.cos((FOUR_PI * (double)x) / denominator));
}
return coefficients;
}
/**
* Creates a 4-term Nutall window
*
* @param length of the window
* @return Nutall window
*/
public static double[] getNutall(int length)
{
double[] coefficients = new double[length];
double denominator = length - 1;
double a0 = 0.355768;
double a1 = 0.487396;
double a2 = 0.144232;
double a3 = 0.012604;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 -
(a1 * Math.cos(TWO_PI * (double)x / denominator)) +
(a2 * Math.cos(FOUR_PI * (double)x / denominator)) -
(a3 * Math.cos(SIX_PI * (double)x / denominator));
}
return coefficients;
}
/**
* Creates a 4-term Blackman-Nutall window
*
* @param length of the window
* @return Blackman-Nutall window
*/
public static double[] getBlackmanNutall(int length)
{
double[] coefficients = new double[length];
double denominator = length - 1;
double a0 = 0.3635819;
double a1 = 0.4891775;
double a2 = 0.1365995;
double a3 = 0.0106411;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 -
(a1 * Math.cos(TWO_PI * (double)x / denominator)) +
(a2 * Math.cos(FOUR_PI * (double)x / denominator)) -
(a3 * Math.cos(SIX_PI * (double)x / denominator));
}
return coefficients;
}
/**
* Creates a 4-term Blackman-Harris window
*
* @param length of the window
* @return Blackman-Harris window
*/
public static double[] getBlackmanHarris4(int length)
{
double[] coefficients = new double[length];
double denominator = length - 1;
double a0 = 0.35875;
double a1 = 0.48829;
double a2 = 0.14128;
double a3 = 0.01168;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 -
(a1 * Math.cos(TWO_PI * (double)x / denominator)) +
(a2 * Math.cos(FOUR_PI * (double)x / denominator)) -
(a3 * Math.cos(SIX_PI * (double)x / denominator));
}
return coefficients;
}
/**
* Creates a 7-term Blackman-Harris window
*
* @param length of the window
* @return Blackman-Harris window
*/
public static double[] getBlackmanHarris7(int length)
{
double[] coefficients = new double[length];
double denominator = length - 1;
double a0 = 0.27105140069342;
double a1 = 0.43329793923448;
double a2 = 0.21812299954311;
double a3 = 0.06592544638803;
double a4 = 0.01081174209837;
double a5 = 0.00077658482522;
double a6 = 0.00001388721735;
return coefficients;
for(int x = 0; x < length; x++)
{
double w = TWO_PI * (double)x / denominator;
coefficients[x] = a0 -
(a1 * Math.cos(w)) +
(a2 * Math.cos(2.0 * w)) -
(a3 * Math.cos(3.0 * w)) +
(a4 * Math.cos(4.0 * w)) -
(a5 * Math.cos(5.0 * w)) +
(a6 * Math.cos(6.0 * w));
}
return coefficients;
}
/**
* Creates a Hamming window
*
* @param length of the window
* @return Hamming window
*/
public static double[] getHamming(int length)
{
double[] coefficients = new double[length];
double a0 = 0.54;
double a1 = 0.46;
double denominator = length;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 - (a1 * Math.cos(TWO_PI * (double)x / denominator));
}
return coefficients;
}
/**
* Creates a Hann window
*
* @param length of the window
* @return Hann window
*/
public static double[] getHann(int length)
{
double[] coefficients = new double[length];
double a0 = 0.5;
double a1 = 0.5;
double denominator = length;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 - (a1 * Math.cos(TWO_PI * (double)x / denominator));
}
return coefficients;
}
/**
* Creates a 4-term Blackman-Harris window
*
* @param length of the window
* @return Blackman-Harris window
*/
public static double[] getFlatTop(int length)
{
double[] coefficients = new double[length];
double denominator = length - 1;
double a0 = 0.215578948;
double a1 = 0.41663158;
double a2 = 0.277263158;
double a3 = 0.083578947;
double a4 = 0.006947368;
for(int x = 0; x < length; x++)
{
coefficients[x] = a0 -
(a1 * Math.cos(TWO_PI * (double)x / denominator)) +
(a2 * Math.cos(FOUR_PI * (double)x / denominator)) -
(a3 * Math.cos(SIX_PI * (double)x / denominator)) +
(a4 * Math.cos(EIGHT_PI * (double)x / denominator));
}
return coefficients;
}
private static double getKaiserBeta(double attenuation)
{
if(attenuation > 50.0)
{
return 0.1102 * (attenuation - 8.7);
}
else if(attenuation >= 21.0)
{
return (0.5842 * Math.pow(attenuation - 21.0, 0.4)) + (0.07886 * (attenuation - 21.0));
}
else
{
return 0.0;
}
}
/**
* Apply the window against an array of float-type samples
*/
public static float[] apply( double[] coefficients, float[] samples )
public static float[] apply(double[] coefficients, float[] samples)
{
for( int x = 0; x < coefficients.length; x++ )
{
samples[ x ] = (float)( samples[ x ] * coefficients[ x ] );
}
return samples;
}
public static float[] apply( WindowType type, float[] samples )
{
double[] coefficients = getWindow( type, samples.length );
return apply( coefficients, samples );
}
/**
* Apply the window against an array of float-type samples
*/
public static double[] apply( double[] coefficients, double[] samples )
{
for( int x = 0; x < coefficients.length; x++ )
{
samples[ x ] = samples[ x ] * coefficients[ x ];
}
return samples;
for(int x = 0; x < coefficients.length; x++)
{
samples[x] = (float)(samples[x] * coefficients[x]);
}
return samples;
}
public static double[] apply( WindowType type, double[] samples )
/**
* Apply the window type against the float array of samples
*
* @param type of window to use
* @param samples to be windowed
* @return windowed samples
*/
public static float[] apply(WindowType type, float[] samples)
{
double[] coefficients = getWindow( type, samples.length );
return apply( coefficients, samples );
double[] coefficients = getWindow(type, samples.length);
return apply(coefficients, samples);
}
/**
* Apply the window against an array of double-type samples
*/
public static double[] apply(double[] coefficients, double[] samples)
{
for(int x = 0; x < coefficients.length; x++)
{
samples[x] = samples[x] * coefficients[x];
}
return samples;
}
/**
* Apply the window type against the double array of samples
*
* @param type of window to use
* @param samples to be windowed
* @return windowed samples
*/
public static double[] apply(WindowType type, double[] samples)
{
double[] coefficients = getWindow(type, samples.length);
return apply(coefficients, samples);
}
/**
* Window types
*/
public enum WindowType
{
BLACKMAN, COSINE, HAMMING, HANNING, NONE;
}
/**
* Utility to log the double arrays
*/
public static String arrayToString( double[] array, boolean breaks )
{
StringBuilder sb = new StringBuilder();
for( int x = 0; x < array.length; x++ )
{
sb.append( x + ":" + array[ x ] );
BLACKMAN("Blackman"),
BLACKMAN_HARRIS_4("Blackman-Harris 4"),
BLACKMAN_HARRIS_7("Blackman-Harris 7"),
BLACKMAN_NUTALL("Blackman-Nutall"),
COSINE("Cosine"),
FLAT_TOP("Flat Top"),
HAMMING("Hamming"),
HANN("Hann"),
NUTALL("Nutall"),
NONE("None");
if( breaks )
{
sb.append( "\n" );
}
else
{
if( x % 8 == 7 )
{
sb.append( "\n" );
}
else
{
sb.append( "\t" );
}
}
}
return sb.toString();
}
private String mLabel;
WindowType(String label)
{
mLabel = label;
}
public String toString()
{
return mLabel;
}
}
}

View File

@@ -0,0 +1,144 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.sample.Broadcaster;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import io.github.dsheirer.source.ISourceEventListener;
import io.github.dsheirer.source.SourceEvent;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.List;
import java.util.concurrent.CopyOnWriteArrayList;
public abstract class AbstractComplexPolyphaseChannelizer implements Listener<ComplexBuffer>, ISourceEventListener
{
private final static Logger mLog = LoggerFactory.getLogger(AbstractComplexPolyphaseChannelizer.class);
private Broadcaster<SourceEvent> mSourceChangeBroadcaster = new Broadcaster();
private List<PolyphaseChannelSource> mChannels = new CopyOnWriteArrayList<>();
private double mSampleRate;
private int mChannelCount;
private double mChannelSampleRate;
/**
* Complex sample polyphase channelizer
*
* @param channelCount
*/
public AbstractComplexPolyphaseChannelizer(double sampleRate, int channelCount)
{
mChannelCount = channelCount;
mSampleRate = sampleRate;
mChannelSampleRate = (double)mSampleRate / (double)mChannelCount;
}
/**
* Input sample rate for this channelizer
* @return sample rate in hertz
*/
public double getSampleRate()
{
return mSampleRate;
}
/**
* Sets the input sample rate for for this channelizer
* @param sampleRate in hertz
*/
public void setSampleRate(double sampleRate)
{
mSampleRate = sampleRate;
mChannelSampleRate = mSampleRate / (double)mChannelCount;
}
/**
* Number of channels being processed by this channelizer.
*/
public int getChannelCount()
{
return mChannelCount;
}
/**
* Output channel sample rate
* @return sample rate in hertz
*/
public double getChannelSampleRate()
{
return mChannelSampleRate;
}
/**
* Dispatches the processed channel samples to any registered polyphase channel outputs.
*
* @param channels containing an array of I/Q samples per channel
*/
protected void dispatch(float[] channels)
{
for(PolyphaseChannelSource channel : mChannels)
{
channel.receiveChannelResults(channels);
}
}
/**
* Adds the polyphase channel source to receive processed output channel samples
*
* @param polyphaseChannelSource
*/
public void addChannel(PolyphaseChannelSource polyphaseChannelSource)
{
if(polyphaseChannelSource != null && !mChannels.contains(polyphaseChannelSource))
{
mChannels.add(polyphaseChannelSource);
mSourceChangeBroadcaster.addListener(polyphaseChannelSource.getSourceEventListener());
}
}
/**
* Removes the polyphase channel source from receiving output channel samples.
*
* @param polyphaseChannelSource
*/
public void removeChannel(PolyphaseChannelSource polyphaseChannelSource)
{
if(polyphaseChannelSource != null && mChannels.contains(polyphaseChannelSource))
{
mChannels.remove(polyphaseChannelSource);
mSourceChangeBroadcaster.removeListener(polyphaseChannelSource.getSourceEventListener());
}
}
/**
* Number of polyphase channels registered to receive sample streams
*/
public int getRegisteredChannelCount()
{
return mChannels.size();
}
@Override
public Listener<SourceEvent> getSourceEventListener()
{
return mSourceChangeBroadcaster;
}
}

View File

@@ -0,0 +1,586 @@
/*******************************************************************************
* *********************************************************************************************************************
* sdr-trunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
* License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with this program.
* If not, see <http://www.gnu.org/licenses/>
* *********************************************************************************************************************
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.source.tuner.channel.TunerChannel;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.List;
public class ChannelCalculator
{
private final static Logger mLog = LoggerFactory.getLogger(ChannelCalculator.class);
/**
* Policy to determine desired corrective action for methods that calculate an index for a frequency or a frequency
* from an index where the frequency may exist either at the overlapping boundary of two indexes, or in the middle
* of the wrap-around (N/2) index.
*/
private enum IndexBoundaryPolicy
{
ADJUST_POSITIVE,
ADJUST_NEGATIVE
}
private double mSampleRate;
private int mChannelCount;
private double mCenterFrequency;
private double mChannelBandwidth;
private double mOversampling = 1.0;
/**
* Calculates channel index(es) from a polyphase channelizer output to source a DDC polyphase channel source
* with a specified channel center frequency and bandwidth.
*
* Polyphase channel indexes are not linear/consecutive owing to the use of the FFT to rotate channels to their
* center frequencies following filtering in the polyphase filter bank. As such, this calculator provides methods
* to determine the min/center/max frequency of each channel that is adjusted for the current center frequency
* value and to determine the channel index(es) that need to be targeted to produce a requested output channel
* with a specific frequency and bandwdith.
*
* @param sampleRate in hertz for the input sample rate from the sample source (ie tuner)
* @param channelCount number of channels in the polyphase filter bank
* @param centerFrequency center frequency for the input sample stream
* @param oversampling set to 1.0 for no oversample, or greater than 1.0 for oversampling
*/
public ChannelCalculator(double sampleRate, int channelCount, double centerFrequency, double oversampling)
{
mSampleRate = sampleRate;
mChannelCount = channelCount;
mCenterFrequency = centerFrequency;
mOversampling = oversampling;
updateChannelBandwidth();;
}
/**
* Calculates the channel bandwidth from the sample rate and channel count.
*/
private void updateChannelBandwidth()
{
mChannelBandwidth = (double)mSampleRate / (double)mChannelCount;
}
/**
* Center frequency for the channelizer
*/
public double getCenterFrequency()
{
return mCenterFrequency;
}
/**
* Updates the center frequency for the channelizer
*/
public void setCenterFrequency(double frequency)
{
mCenterFrequency = frequency;
}
/**
* Number of channels being channelized by the polyphase channelizer
*/
public int getChannelCount()
{
return mChannelCount;
}
/**
* Updates the number of channels being channelized by the polyphase channelizer
*/
public void setChannelCount(int channelCount)
{
mChannelCount = channelCount;
updateChannelBandwidth();
}
/**
* Identifies the index of the wrap-around channel meaning that the lower half of the channel exists at the upper
* end of the spectrum and the upper half of the channel exists at the lower end of the spectrum
*
* @return wrap-around channel index
*/
public int getWrapAroundIndex()
{
return mChannelCount / 2;
}
/**
* PolyphaseChannelManager input sample rate
* @return input sample rate in hertz
*/
public double getSampleRate()
{
return mSampleRate;
}
/**
* Updates the channel count using the sample rate argument to calculate the number of channels from the
* already established channel bandwidth.
*
* @param sampleRate to use for channel count
*/
public void setSampleRate(double sampleRate)
{
setChannelCount((int)(sampleRate / getChannelBandwidth()));
}
/**
* Channel bandwidth. Note: this value may be smaller than the channel sample rate when oversampling
* is being used.
*/
public double getChannelBandwidth()
{
return mChannelBandwidth;
}
/**
* Half of the channel bandwidth.
*/
public double getHalfChannelBandwidth()
{
return mChannelBandwidth / 2.0;
}
/**
* Oversampling rate produced by the channelizer.
*
* @return 1.0 for no oversampling, or a value greater than 1.0 when oversampling is applied.
*/
public double getOversampling()
{
return mOversampling;
}
/**
* Sets the oversampling rate being applied by the polyphase channelizer
* @param oversampling rate where 1.0 indicates no oversampling, and a value greater than 1.0 indicates
* that the polyphase channelizer is oversampling each channel.
*/
public void setOversampling(double oversampling)
{
mOversampling = oversampling;
}
/**
* Channel sample rate. This value is the product of the filtered channel bandwidth times the
* oversampling rate.
*/
public double getChannelSampleRate()
{
return (double)getChannelBandwidth() * getOversampling();
}
/**
* PolyphaseChannelManager input minimum frequency
* @return minimum frequency in hertz
*/
public long getMinimumFrequency()
{
return (long)(getCenterFrequency() - (getSampleRate() / 2.0));
}
/**
* PolyphaseChannelManager input maximum frequency
* @return maximum frequency in hertz
*/
public long getMaximumFrequency()
{
return (long)(getCenterFrequency() + (getSampleRate() / 2.0));
}
/**
* Identifies the channelizer channel index(es) required to produce a tuner channel with the specified center
* frequency and channel bandwidth.
*
* @param tunerChannel specifying the desired channel frequency and bandwidth
*
* @return an array of one or more polyphase channelizer index(es) to target in order to produce an output channel
* for the desired frequency and bandwidth.
*
* @throws IllegalArgumentException if the tuner channel's minimum or maximum frequency is not contained
* within the channelizer's minimum or maximum frequency range
*/
public List<Integer> getChannelIndexes(TunerChannel tunerChannel) throws IllegalArgumentException
{
if(tunerChannel.getMinFrequency() < getMinimumFrequency() ||
tunerChannel.getMaxFrequency() > getMaximumFrequency())
{
throw new IllegalArgumentException("Requested channel cannot be provided by this " +
"channelizer. Requested channel [" + tunerChannel.getMinFrequency() + " - " +
tunerChannel.getMaxFrequency() + "] exceeds current channelizer frequency range [" +
getMinimumFrequency() + " - " + getMaximumFrequency() + "]");
}
int minIndex = getIndexForFrequency(tunerChannel.getMinFrequency(), IndexBoundaryPolicy.ADJUST_POSITIVE);
int maxIndex = getIndexForFrequency(tunerChannel.getMaxFrequency(), IndexBoundaryPolicy.ADJUST_NEGATIVE);
//Detect the unsourceable situation where the wrap around index would be included in both the minimum index
//and the maximum index, meaning that the tuner channel bandwidth is nearly the same as the input sample rate.
if(minIndex == getWrapAroundIndex() && maxIndex == getWrapAroundIndex())
{
throw new IllegalArgumentException("Requested tuner channel cannot be provided. Requested bandwidth is " +
"within two channel bandwidths of the sample rate.");
}
List<Integer> indexes = new ArrayList<>();
indexes.add(minIndex);
if(minIndex != maxIndex)
{
//Ensure our min/max indexes are within the correct bounds 0 <> N-1 to safeguard the while loop below
if(minIndex < 0 || minIndex >= getChannelCount() || maxIndex < 0 || maxIndex >= getChannelCount())
{
throw new IllegalArgumentException("Something went wrong while calculating the min and max polyphase " +
"channel indexes for tuner channel: " + tunerChannel);
}
int pointer = minIndex + 1;
if(pointer >= getChannelCount())
{
pointer -= getChannelCount();
}
while(pointer != maxIndex)
{
indexes.add(pointer);
pointer++;
if(pointer >= getChannelCount())
{
pointer -= getChannelCount();
}
}
indexes.add(maxIndex);
}
return indexes;
}
/**
* Calculates the index for the specified frequency. Certain frequencies can exist as both the maximum frequency
* of one channel and the minimum frequency of an adjacent channel. The index boundary policy identifies your
* choice to select the higher frequency channel index (ADJUST_POSITIVE) or the lower frequency channel
* index (ADJUST_NEGATIVE).
*
* @param frequency to determine channel index for
* @param indexBoundaryPolicy to determine the behavior when the frequency overlaps two channel boundaries
* @return index for the frequency
*/
public int getIndexForFrequency(long frequency, IndexBoundaryPolicy indexBoundaryPolicy)
{
double offset = frequency - getCenterFrequency();
if(Math.abs(offset) < getHalfChannelBandwidth())
{
return 0; //Center index
}
if(offset > 0)
{
offset += getHalfChannelBandwidth();
}
else
{
offset -= getHalfChannelBandwidth();
}
int indexOffset = (int)(offset / getChannelBandwidth());
if(indexOffset < 0)
{
indexOffset += getChannelCount();
}
if(indexBoundaryPolicy == IndexBoundaryPolicy.ADJUST_POSITIVE &&
isOverlapFrequency(frequency, indexOffset, indexOffset + 1))
{
indexOffset = normalize(indexOffset + 1);
}
else if(indexBoundaryPolicy == IndexBoundaryPolicy.ADJUST_NEGATIVE &&
isOverlapFrequency(frequency, indexOffset - 1, indexOffset))
{
indexOffset = normalize(indexOffset - 1);
}
return indexOffset;
}
/**
* Indicates if the frequency is an overlapping frequency between channel index 1 and channel index 2 where the
* maximum frequency of index 1 is the same as the minimum frequency of index 2 and both of these frequencies are
* the same as the frequency argument.
*
* Note: both index 1 and index 2 will be normalized to the range of 0 <> N-1 prior to determining overlap
*
* @param frequency to test for overlap between index 1 and 2
* @param index1 which is the lower channel index located next to index 2.
* @param index2 which is the higher channel index located next to index 1.
* @return true if frequency = index1 maximum frequency = index2 minimum frequency
*/
public boolean isOverlapFrequency(long frequency, int index1, int index2)
{
index1 = normalize(index1);
index2 = normalize(index2);
//Check that these two indexes are side by side
int delta = normalize(index2 - index1);
if(delta != 1)
{
return false;
}
long index1MaximumFrequency = (long)getIndexMaximumFrequency(index1, IndexBoundaryPolicy.ADJUST_POSITIVE);
long index2MinimumFrequency = (long)getIndexMinimumFrequency(index2, IndexBoundaryPolicy.ADJUST_NEGATIVE);
if(index1 == getWrapAroundIndex())
{
index1MaximumFrequency = (long)getIndexMaximumFrequency(index1, IndexBoundaryPolicy.ADJUST_NEGATIVE);
}
if(index2 == getWrapAroundIndex())
{
index2MinimumFrequency = (long)getIndexMinimumFrequency(index2, IndexBoundaryPolicy.ADJUST_POSITIVE);
}
return frequency == index1MaximumFrequency && frequency == index2MinimumFrequency;
}
/**
* Normalizes the index to the range of 0 <> N-1 channel indexes</>
*/
private int normalize(int index)
{
while(index < 0)
{
index += getChannelCount();
}
while(index >= getChannelCount())
{
index -= getChannelCount();
}
return index;
}
/**
* Calculates the center frequency for the channel identified by the index in the channel results array.
*
* @param index to determine the channel center frequency
* @return center frequency for the indexed channel
*/
public double getIndexCenterFrequency(int index, IndexBoundaryPolicy indexBoundaryPolicy)
{
if(index < 0 || index >= getChannelCount())
{
throw new IllegalArgumentException("Illegal channel index");
}
int wrapAroundIndex = getWrapAroundIndex();
if(index == wrapAroundIndex)
{
if(indexBoundaryPolicy == IndexBoundaryPolicy.ADJUST_POSITIVE)
{
return getCenterFrequency() + (index * getChannelBandwidth());
}
else
{
return getCenterFrequency() - (index * getChannelBandwidth());
}
}
else if(index < wrapAroundIndex)
{
return getCenterFrequency() + (index * getChannelBandwidth());
}
else
{
return getCenterFrequency() - ((getChannelCount() - index) * getChannelBandwidth());
}
}
/**
* Calculates the minimum frequency band edge for the channel identified by the index in the channel results
* array.
* @param index to calculate the minimum band edge frequency for the channel.
* @param indexBoundaryPolicy for the wrap around index to indicate if this method should return the minimum
* frequency of the upper spectrum half of the wrap around channel (ADJUST_POSITIVE) or the minimum frequency of the
* lower spectrum half of the wrap around channel (ADJUST_NEGATIVE)
* @return frequency of the minimum band edge
*/
public double getIndexMinimumFrequency(int index, IndexBoundaryPolicy indexBoundaryPolicy)
{
if(index < 0 || index >= getChannelCount())
{
throw new IllegalArgumentException("Can't determine minimum frequency - illegal channel index [" +
index + "]");
}
int wrapAroundIndex = getWrapAroundIndex();
if(index == wrapAroundIndex)
{
if(indexBoundaryPolicy == IndexBoundaryPolicy.ADJUST_POSITIVE)
{
return getCenterFrequency() + ((double)index * getChannelBandwidth()) - getHalfChannelBandwidth();
}
else
{
return getIndexCenterFrequency(index, indexBoundaryPolicy);
}
}
else if(index <= wrapAroundIndex)
{
return getCenterFrequency() + ((double)index * getChannelBandwidth()) - getHalfChannelBandwidth();
}
else
{
return getCenterFrequency() - ((double)(getChannelCount() - index) * getChannelBandwidth()) -
getHalfChannelBandwidth();
}
}
/**
* Calculates the maximum frequency band edge for the channel identified by the index in the channel results
* array.
* @param index to calculate the maximum band edge frequency for the channel.
* @param indexBoundaryPolicy for the wrap around index to indicate if this method should return the maximum
* frequency of the upper spectrum half of the wrap around channel (ADJUST_POSITIVE) or the maximum frequency of the
* lower spectrum half of the wrap around channel (ADJUST_NEGATIVE)
* @return frequency of the maximum band edge
*/
public double getIndexMaximumFrequency(int index, IndexBoundaryPolicy indexBoundaryPolicy)
{
if(index < 0 || index >= getChannelCount())
{
throw new IllegalArgumentException("Can't determine maximum frequency - illegal channel index [" +
index + "]");
}
int wrapAroundIndex = getWrapAroundIndex();
if(index == wrapAroundIndex)
{
if(indexBoundaryPolicy == IndexBoundaryPolicy.ADJUST_POSITIVE)
{
return getIndexCenterFrequency(index, indexBoundaryPolicy);
}
else
{
return getCenterFrequency() - ((double)index * getChannelBandwidth() - getHalfChannelBandwidth());
}
}
else if(index <= wrapAroundIndex)
{
return getCenterFrequency() + ((double)index * getChannelBandwidth()) + getHalfChannelBandwidth();
}
else
{
return getCenterFrequency() - ((double)(getChannelCount() - index) * getChannelBandwidth()) +
getHalfChannelBandwidth();
}
}
/**
* Calculates the center frequency for a channel that will be synthesized from the channel indexes.
* @param indexes that are contiguous in the spectrum
* @return center frequency in hertz
*/
public long getCenterFrequencyForIndexes(List<Integer> indexes)
{
if(indexes.isEmpty())
{
throw new IllegalArgumentException("Indexes cannot be empty");
}
int middleIndex = indexes.get(indexes.size() / 2);
if(indexes.size() % 2 == 0)
{
//Even number of indexes
return (long)getIndexMinimumFrequency(middleIndex, IndexBoundaryPolicy.ADJUST_POSITIVE);
}
else
{
//Odd number of indexes
return (long)getIndexCenterFrequency(middleIndex, IndexBoundaryPolicy.ADJUST_POSITIVE);
}
}
/**
* String representation of the channel calculator settings
*/
public String toString()
{
StringBuilder sb = new StringBuilder();
sb.append("Channels:").append(getChannelCount());
sb.append(" Sample Rate:").append(getSampleRate());
sb.append(" Channel Bandwidth:").append(getChannelBandwidth());
sb.append(" Channel Rate:").append(getChannelSampleRate());
sb.append(" Min:").append(getMinimumFrequency());
sb.append(" Max:").append(getMaximumFrequency());
return sb.toString();
}
public static void main(String[] args)
{
mLog.debug("Starting ....");
int channelBandwidth = 12500;
int sampleRate = 100000; //Airspy
int channelCount = sampleRate / channelBandwidth;
double centerFrequency = 455000000;
ChannelCalculator calculator = new ChannelCalculator(sampleRate, channelCount, centerFrequency, 2.0);
mLog.debug(calculator.toString());
DecimalFormat decimalFormat = new DecimalFormat("000,000,000.0");
int wrapAroundIndex = calculator.getWrapAroundIndex();
IndexBoundaryPolicy indexBoundaryPolicy = IndexBoundaryPolicy.ADJUST_NEGATIVE;
for(int index = wrapAroundIndex; index < calculator.getChannelCount(); index++)
{
mLog.debug("Index:" + index + " Min:" + decimalFormat.format(calculator.getIndexMinimumFrequency(index, indexBoundaryPolicy)) +
" Center:" + decimalFormat.format(calculator.getIndexCenterFrequency(index, indexBoundaryPolicy)) +
" Max:" + decimalFormat.format(calculator.getIndexMaximumFrequency(index, indexBoundaryPolicy)));
}
indexBoundaryPolicy = IndexBoundaryPolicy.ADJUST_POSITIVE;
for(int index = 0; index <= wrapAroundIndex; index++)
{
mLog.debug("Index:" + index + " Min:" + decimalFormat.format(calculator.getIndexMinimumFrequency(index, indexBoundaryPolicy)) +
" Center:" + decimalFormat.format(calculator.getIndexCenterFrequency(index, indexBoundaryPolicy)) +
" Max:" + decimalFormat.format(calculator.getIndexMaximumFrequency(index, indexBoundaryPolicy)));
}
mLog.debug("Finished");
}
}

View File

@@ -0,0 +1,344 @@
/*******************************************************************************
* *********************************************************************************************************************
* sdr-trunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
* License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with this program.
* If not, see <http://www.gnu.org/licenses/>
* *********************************************************************************************************************
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.SampleType;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import io.github.dsheirer.settings.SettingsManager;
import io.github.dsheirer.source.ISourceEventProcessor;
import io.github.dsheirer.source.SourceEvent;
import io.github.dsheirer.source.SourceException;
import io.github.dsheirer.source.tuner.channel.TunerChannel;
import io.github.dsheirer.source.tuner.channel.TunerChannelSource;
import io.github.dsheirer.source.tuner.configuration.TunerConfigurationModel;
import io.github.dsheirer.source.tuner.test.TestTuner;
import io.github.dsheirer.spectrum.DFTProcessor;
import io.github.dsheirer.spectrum.DFTSize;
import io.github.dsheirer.spectrum.SpectrumPanel;
import io.github.dsheirer.spectrum.converter.ComplexDecibelConverter;
import net.miginfocom.swing.MigLayout;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import javax.swing.JFrame;
import javax.swing.JLabel;
import javax.swing.JPanel;
import javax.swing.JSpinner;
import javax.swing.SpinnerNumberModel;
import javax.swing.event.ChangeEvent;
import javax.swing.event.ChangeListener;
import java.awt.Dimension;
import java.awt.EventQueue;
public class ChannelizerViewer extends JFrame
{
private final static Logger mLog = LoggerFactory.getLogger(ChannelizerViewer.class);
private static final int CHANNEL_BANDWIDTH = 12500;
private static final int CHANNEL_FFT_FRAME_RATE = 20;
private SettingsManager mSettingsManager = new SettingsManager(new TunerConfigurationModel());
private JPanel mPrimaryPanel;
private JPanel mControlPanel;
private PrimarySpectrumPanel mPrimarySpectrumPanel;
private ChannelArrayPanel mChannelPanel;
private int mChannelCount;
private int mChannelsPerRow;
private long mToneFrequency;
private DFTSize mMainPanelDFTSize = DFTSize.FFT32768;
private DFTSize mChannelPanelDFTSize = DFTSize.FFT01024;
private TestTuner mTestTuner;
/**
* GUI Test utility for researching polyphase channelizers.
*
* @param channelsPerRow
*/
public ChannelizerViewer(int channelsPerRow)
{
mTestTuner = new TestTuner();
mChannelCount = (int)(mTestTuner.getTunerController().getSampleRate() / CHANNEL_BANDWIDTH);
mChannelsPerRow = channelsPerRow;
mToneFrequency = mTestTuner.getTunerController().getFrequency() + 19200;
mTestTuner.getTunerController().setToneFrequency(mToneFrequency);
init();
}
private void init()
{
setTitle("TunerChannelizer Viewer");
setSize(1000, 800);
setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
setLayout(new MigLayout("insets 0 0 0 0", "[grow,fill]", "[grow,fill]"));
setLocationRelativeTo(null);
add(getPrimaryPanel());
}
private JPanel getPrimaryPanel()
{
if(mPrimaryPanel == null)
{
mPrimaryPanel = new JPanel();
mPrimaryPanel.setLayout(new MigLayout("insets 0 0 0 0", "[grow,fill]", "[][][grow,fill]"));
mPrimaryPanel.add(getSpectrumPanel(), "wrap");
mPrimaryPanel.add(getControlPanel(), "wrap");
mPrimaryPanel.add(getChannelArrayPanel());
}
return mPrimaryPanel;
}
private PrimarySpectrumPanel getSpectrumPanel()
{
if(mPrimarySpectrumPanel == null)
{
mPrimarySpectrumPanel = new PrimarySpectrumPanel(mSettingsManager,
mTestTuner.getTunerController().getSampleRate());
mPrimarySpectrumPanel.setPreferredSize(new Dimension(1000, 200));
mPrimarySpectrumPanel.setDFTSize(mMainPanelDFTSize);
mTestTuner.getSourceManager().addComplexBufferListener(mPrimarySpectrumPanel);
}
return mPrimarySpectrumPanel;
}
private JPanel getControlPanel()
{
if(mControlPanel == null)
{
mControlPanel = new JPanel();
mControlPanel.setLayout(new MigLayout("insets 0 0 0 0", "", ""));
mControlPanel.add(new JLabel("Tone:"), "align left");
int maximumFrequencyOffset = (int)(mTestTuner.getTunerController().getSampleRate() / 2);
long minimumFrequency = mTestTuner.getTunerController().getFrequency() - maximumFrequencyOffset;
long maximumFrequency = mTestTuner.getTunerController().getFrequency() + maximumFrequencyOffset;
SpinnerNumberModel model = new SpinnerNumberModel(mToneFrequency, minimumFrequency, maximumFrequency,
100 );
model.addChangeListener(new ChangeListener()
{
@Override
public void stateChanged(ChangeEvent e)
{
int toneFrequency = model.getNumber().intValue();
mTestTuner.getTunerController().setToneFrequency(toneFrequency);
}
});
JSpinner spinner = new JSpinner(model);
mControlPanel.add(spinner);
mControlPanel.add(new JLabel("Hz"));
mControlPanel.add(new JLabel("Sample Rate: " + mTestTuner.getTunerController().getSampleRate() + " Hz"), "push,align right");
mControlPanel.add(new JLabel("Max: " + (mTestTuner.getTunerController().getSampleRate() / 2) + " Hz"), "push,align right");
mControlPanel.add(new JLabel("Channels: " + mChannelCount), "push,align right");
}
return mControlPanel;
}
private ChannelArrayPanel getChannelArrayPanel()
{
if(mChannelPanel == null)
{
mChannelPanel = new ChannelArrayPanel();
}
return mChannelPanel;
}
public class ChannelArrayPanel extends JPanel
{
public ChannelArrayPanel()
{
int bufferSize = CHANNEL_BANDWIDTH / CHANNEL_FFT_FRAME_RATE;
if(bufferSize % 2 == 1)
{
bufferSize++;
}
init();
}
private void init()
{
setLayout(new MigLayout("insets 0 0 0 0", "fill", "fill"));
double spectralBandwidth = mTestTuner.getTunerController().getSampleRate();
double halfSpectralBandwidth = spectralBandwidth / 2.0;
int channelToLog = -1;
for(int x = 0; x < mChannelCount; x++)
{
//place the channels left and right of 10.0 MHz
long frequency = (long)(mTestTuner.getTunerController().getFrequency() +
((x * CHANNEL_BANDWIDTH) - halfSpectralBandwidth) + (CHANNEL_BANDWIDTH / 2));
mLog.debug("Channel " + x + " Frequency: " + frequency);
ChannelPanel channelPanel = new ChannelPanel(mSettingsManager, CHANNEL_BANDWIDTH * 2, frequency, CHANNEL_BANDWIDTH, (x == channelToLog));
channelPanel.setDFTSize(mChannelPanelDFTSize);
if(x % mChannelsPerRow == mChannelsPerRow - 1)
{
add(channelPanel, "grow,push,wrap 2px");
}
else
{
add(channelPanel, "grow,push");
}
}
}
}
public class PrimarySpectrumPanel extends JPanel implements Listener<ComplexBuffer>, ISourceEventProcessor
{
private DFTProcessor mDFTProcessor = new DFTProcessor(SampleType.COMPLEX);
private ComplexDecibelConverter mComplexDecibelConverter = new ComplexDecibelConverter();
private SpectrumPanel mSpectrumPanel;
public PrimarySpectrumPanel(SettingsManager settingsManager, double sampleRate)
{
setLayout(new MigLayout("insets 0 0 0 0", "[grow,fill]", "[grow,fill]"));
mSpectrumPanel = new SpectrumPanel(settingsManager);
add(mSpectrumPanel);
mDFTProcessor.addConverter(mComplexDecibelConverter);
mDFTProcessor.process(SourceEvent.sampleRateChange(sampleRate));
mComplexDecibelConverter.addListener(mSpectrumPanel);
}
public void setDFTSize(DFTSize dftSize)
{
mDFTProcessor.setDFTSize(dftSize);
}
@Override
public void receive(ComplexBuffer complexBuffer)
{
mDFTProcessor.receive(complexBuffer);
}
@Override
public void process(SourceEvent event) throws SourceException
{
mLog.debug("Source Event! Add handler support for this to channelizer viewer");
}
}
public class ChannelPanel extends JPanel implements Listener<ComplexBuffer>, ISourceEventProcessor
{
private TunerChannelSource mSource;
private DFTProcessor mDFTProcessor = new DFTProcessor(SampleType.COMPLEX);
private ComplexDecibelConverter mComplexDecibelConverter = new ComplexDecibelConverter();
private SpectrumPanel mSpectrumPanel;
private long mFrequency;
public ChannelPanel(SettingsManager settingsManager, double sampleRate, long frequency, int bandwidth, boolean enableLogging)
{
setLayout(new MigLayout("insets 0 0 0 0", "[grow,fill]", "[grow,fill]"));
mSpectrumPanel = new SpectrumPanel(settingsManager);
add(mSpectrumPanel);
mDFTProcessor.addConverter(mComplexDecibelConverter);
mDFTProcessor.process(SourceEvent.sampleRateChange(sampleRate));
mComplexDecibelConverter.addListener(mSpectrumPanel);
mFrequency = frequency;
TunerChannel tunerChannel = new TunerChannel(frequency, bandwidth);
mSource = mTestTuner.getSourceManager().getSource(tunerChannel);
if(mSource != null)
{
if(enableLogging)
{
mSource.setListener(new Listener<ComplexBuffer>()
{
@Override
public void receive(ComplexBuffer complexBuffer)
{
mLog.debug("Sending samples to DFT processor: " + complexBuffer.getSamples().length);
mDFTProcessor.receive(complexBuffer);
}
});
}
else
{
mSource.setListener(mDFTProcessor);
}
mSource.start();
}
else
{
mLog.error("Couldn't get a source from the tuner for frequency: " + frequency);
}
}
public TunerChannelSource getSource()
{
return mSource;
}
public void setDFTSize(DFTSize dftSize)
{
mDFTProcessor.setDFTSize(dftSize);
}
@Override
public void receive(ComplexBuffer complexBuffer)
{
mLog.debug("Got a buffer for the channel dft processor");
mDFTProcessor.receive(complexBuffer);
}
@Override
public void process(SourceEvent event) throws SourceException
{
mLog.debug("Source Event! Add handler support for this to channelizer viewer");
}
}
public static void main(String[] args)
{
int channelsPerRow = 4;
final ChannelizerViewer frame = new ChannelizerViewer(channelsPerRow);
EventQueue.invokeLater(new Runnable()
{
@Override
public void run()
{
frame.setVisible(true);
}
});
}
}

View File

@@ -0,0 +1,171 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.dsp.filter.FilterFactory;
import io.github.dsheirer.dsp.filter.Window;
import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter;
import io.github.dsheirer.dsp.mixer.Oscillator;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class ComplexPolyphaseChannelizer extends AbstractComplexPolyphaseChannelizer
{
private final static Logger mLog = LoggerFactory.getLogger(ComplexPolyphaseChannelizer.class);
private RealFIRFilter[] mInphaseFilters;
private RealFIRFilter[] mQuadratureFilters;
private int mFilterPointer = 0;
private float[] mFilteredSamples;
private FloatFFT_1D mFFT;
/**
* Maximally Decimated Polyphase Filter Bank (MDPFB) channelizer that divides the input frequency band into
* equal bandwidth channels.
*
* @param coefficients of a prototype low-pass filter designed for the inbound sample rate with a cutoff frequency
* equal to the channel bandwidth (sample rate / filters). If you need to synthesize (combine two or more
* filter outputs) a new bandwidth signal from the outputs of this filter, then the filter should be designed
* as a nyquist filter with -6 dB attenuation at the channel bandwidth cutoff frequency
*
* @param channelCount - number of channels to output.
* @param channelSampleRate - channel sample rate (after decimation)
*/
public ComplexPolyphaseChannelizer(float[] coefficients, int channelCount, int channelSampleRate)
{
super(channelCount, channelSampleRate);
initFilters(coefficients);
}
@Override
public void receive(ComplexBuffer complexBuffer)
{
float[] samples = complexBuffer.getSamples();
for(int x = 0; x < samples.length; x += 2)
{
filter(samples[x], samples[x + 1]);
}
}
protected void filter(float inphase, float quadrature)
{
int index = (getChannelCount() - mFilterPointer - 1) * 2;
mFilteredSamples[index] = mInphaseFilters[mFilterPointer].filter(inphase);
mFilteredSamples[index + 1] = mQuadratureFilters[mFilterPointer].filter(quadrature);
if(mFilterPointer == 0)
{
calculate();
mFilterPointer = getChannelCount() - 1;
}
else
{
mFilterPointer--;
}
}
private void calculate()
{
float[] samples = new float[getChannelCount() * 2];
System.arraycopy(mFilteredSamples, 0, samples, 0, mFilteredSamples.length);
//FFT is executed in-place where the output overwrites the input
mFFT.complexForward(samples);
dispatch(samples);
}
/**
* Distributes the filter taps to each of the polyphase filters
*
* @param taps
*/
private void initFilters(float[] taps)
{
int tapCount = (int)Math.ceil((double)taps.length / (double)getChannelCount());
float[][] coefficients = new float[getChannelCount()][tapCount];
for(int tap = 0; tap < tapCount; tap++)
{
for(int filter = 0; filter < getChannelCount(); filter++)
{
int index = tap * getChannelCount() + filter;
if(index < taps.length)
{
coefficients[filter][tap] = taps[index];
}
}
}
mInphaseFilters = new RealFIRFilter[getChannelCount()];
mQuadratureFilters = new RealFIRFilter[getChannelCount()];
for(int filter = 0; filter < getChannelCount(); filter++)
{
mInphaseFilters[filter] = new RealFIRFilter(coefficients[filter], 1.0f);
mQuadratureFilters[filter] = new RealFIRFilter(coefficients[filter], 1.0f);
}
mFilteredSamples = new float[getChannelCount() * 2];
mFFT = new FloatFFT_1D(getChannelCount());
mFilterPointer = getChannelCount() -1;
}
public static void main(String[] args)
{
mLog.debug("Starting ...");
int symbolRate = 4800;
int samplesPerSymbol = 2;
int channels = 800;
int channelBandwidth = 12500;
int sampleRate = channels * channelBandwidth;
float[] taps = FilterFactory.getLowPass( 10000000, 12500, 5600, Window.WindowType.BLACKMAN_HARRIS_7 );
int channelCenter = (int)(12500 * 1.5);
StringBuilder sb = new StringBuilder();
sb.append("\nPolyphase TunerChannelizer\n");
sb.append("Sample Rate:" + sampleRate + " Channels:" + channels + " Channel Rate:" + channelBandwidth + "\n");
sb.append("Tap Count:" + taps.length + "\n");
sb.append("Channel:" + channelCenter);
mLog.debug(sb.toString());
ComplexPolyphaseChannelizer channelizer = new ComplexPolyphaseChannelizer(taps, 8, channelBandwidth);
Oscillator oscillator = new Oscillator(channelCenter, 2 * sampleRate);
for(int x = 0; x < 2000; x++)
{
channelizer.filter(oscillator.getComplex().inphase(), oscillator.getComplex().quadrature());
oscillator.rotate();
}
mLog.debug("Finished!");
}
}

View File

@@ -0,0 +1,349 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.dsp.filter.FilterFactory;
import io.github.dsheirer.dsp.filter.Window.WindowType;
import io.github.dsheirer.dsp.filter.design.FilterDesignException;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChannelizer
{
private final static Logger mLog = LoggerFactory.getLogger(ComplexPolyphaseChannelizerM2.class);
private static final int DEFAULT_MINIMUM_CHANNEL_BANDWIDTH = 12500;
private float[][] mSamples;
private float[][] mCoefficients;
private int[][] mSampleIndexMap;
private int mColumnPointer;
private int mRowPointer;
private int mBlockSize;
private FloatFFT_1D mFFT;
/**
* Non-Maximally Decimated Polyphase Filter Bank (NMDPFB) channelizer that divides the input frequency band into
* equal bandwidth channels that are each oversampled by 2x for output.
*
* @param taps of a low-pass filter designed for the inbound sample rate with a cutoff frequency
* equal to the channel bandwidth (sample rate / filters). If you need to synthesize (combine two or more
* filter outputs) a new bandwidth signal from the outputs of this filter, then the filter should be designed
* as a nyquist filter with -6 dB attenuation at the channel bandwidth cutoff frequency
*
* @param sampleRate of the incoming sample stream
* @param channelCount - number of filters/channels to output. Since this filter bank oversamples each filter
* output, this number must be even (divisible by the 2x oversample rate).
*/
public ComplexPolyphaseChannelizerM2(float[] taps, int sampleRate, int channelCount)
{
super(sampleRate, channelCount);
if(channelCount % 2 != 0)
{
throw new IllegalArgumentException("Channel count must be an even multiple of the oversample rate (2)");
}
init(taps);
}
/**
* Non-Maximally Decimated Polyphase Filter Bank (NMDPFB) channelizer that divides the input frequency band into
* equal bandwidth channels that are each oversampled by 2x for output.
*
* Note: this constructor uses a default minimum channel bandwidth of 12.5 kHz. The number of channels is calculate
* from the sample rate and this default minimum channel bandwidth. The resulting per-channel bandwidth will be
* between 12.5 and 25.0 kHz with the output per-channel sample rate of 2x the bandwidth.
*
* @param sampleRate to be channelized.
* @param tapsPerFilter to use when designing the filter
*/
public ComplexPolyphaseChannelizerM2(double sampleRate, int tapsPerFilter) throws FilterDesignException
{
super(sampleRate, (int)(sampleRate / DEFAULT_MINIMUM_CHANNEL_BANDWIDTH));
float[] filterTaps = FilterFactory.getSincChannelizer(getChannelSampleRate(), getChannelCount(),
tapsPerFilter, WindowType.BLACKMAN_HARRIS_7, true);
init(filterTaps);
}
/**
* Updates this channelizer to use the new sample rate. This method creates a new filter suitable for the
* sample rate and reinitializes all internal data structures to prepare for processing the new sample rate.
* @param sampleRate in hertz
*/
@Override
public void setSampleRate(double sampleRate)
{
try
{
super.setSampleRate(sampleRate);
float[] filterTaps = FilterFactory.getSincChannelizer(getChannelSampleRate(), getChannelCount(),
15, WindowType.BLACKMAN_HARRIS_7, true);
init(filterTaps);
}
catch(FilterDesignException fde)
{
throw new IllegalArgumentException("Cannot create a channelizer filter for the specified sample rate [" +
sampleRate + "]");
}
}
/**
* Receives the complex sample buffer and processes the results through the channelizer.
*/
@Override
public void receive(ComplexBuffer complexBuffer)
{
float[] samples = complexBuffer.getSamples();
int samplesPointer = 0;
while(samplesPointer < samples.length)
{
int lengthToCopy = mBlockSize - (mRowPointer % mBlockSize);
if((samplesPointer + lengthToCopy) > samples.length)
{
lengthToCopy = samples.length - samplesPointer;
}
System.arraycopy(samples, samplesPointer, mSamples[mColumnPointer], mRowPointer, lengthToCopy);
samplesPointer += lengthToCopy;
mRowPointer += lengthToCopy;
if(mRowPointer % mBlockSize == 0)
{
process();
}
}
}
/**
* Processes the sample buffer for each new block of sample data that is loaded and distributes the results to any
* registered channel listeners.
*/
private void process()
{
int sampleRowCount = getChannelCount() * 2;
float[] processed = new float[sampleRowCount];
//Top of block processing
if(mRowPointer >= sampleRowCount)
{
int[] sampleIndexMap = mSampleIndexMap[mColumnPointer];
for(int channel = 0; channel < getChannelCount(); channel++)
{
int outputIndex = 2 * (getChannelCount() - channel - 1);
int sampleIndex = 2 * channel;
float[] filter = mCoefficients[channel];
processed[outputIndex] = filter((sampleIndex), mCoefficients[channel], sampleIndexMap);
processed[outputIndex + 1] = filter((sampleIndex + 1), mCoefficients[channel], sampleIndexMap);
}
mRowPointer = 0;
mColumnPointer--;
if(mColumnPointer < 0)
{
mColumnPointer = mSamples.length - 1;
}
}
//Mid-point processing
else
{
int[] sampleIndexMap = mSampleIndexMap[mColumnPointer];
int half = getChannelCount() / 2;
for(int channel = 0; channel < half; channel++)
{
int outputIndex = 2 * (getChannelCount() - channel - 1);
int sampleIndex = 2 * channel;
int filterIndex = channel + half;
processed[outputIndex] = filter(sampleIndex, mCoefficients[filterIndex], sampleIndexMap);
processed[outputIndex + 1] = filter(sampleIndex + 1, mCoefficients[filterIndex], sampleIndexMap);
}
//The second half of the sample buffer was loaded using the previous column pointer - adjust for that now
sampleIndexMap = mSampleIndexMap[(mColumnPointer == (mSamples.length - 1) ? 0 : mColumnPointer + 1)];
for(int channel = half; channel < getChannelCount(); channel++)
{
int outputIndex = 2 * (getChannelCount() - channel - 1);
int sampleIndex = 2 * channel;
int filterIndex = channel - half;
processed[outputIndex] = filter(sampleIndex, mCoefficients[filterIndex], sampleIndexMap);
processed[outputIndex + 1] = filter(sampleIndex + 1, mCoefficients[filterIndex], sampleIndexMap);
}
}
//Rotate each of the channels to the correct phase using the IFFT
mFFT.complexInverse(processed, true);
dispatch(processed);
}
/**
* Calculates the filtered output for polyphase filter arm.
*
* @param sampleRow index pointing to a row in the mSamples array to filter
* @param coefficients to use in filtering the samples
* @param sampleIndexMap to translate the column index offset based on the current sample column pointer
* @return filtered sample for the row
*/
private float filter(int sampleRow, float[] coefficients, int[] sampleIndexMap)
{
float accumulator = 0.0f;
for(int column = 0; column < mSamples.length; column++)
{
accumulator += mSamples[sampleIndexMap[column]][sampleRow] * coefficients[column];
}
return accumulator;
}
/**
* Initializes the channelizer structures. Distributes the prototype filter coefficients to each of the polyphase
* filters, sets up the sample buffer array and pointers and initializes the FFT.
*
* These structures resemble those described by Fred Harris in Multirate Signal Processing for Communications
* Systems, p230-233. However, instead of swapping the input sample buffers and swapping the output sample array
* on every other processing block, we simply swap the filter coefficient banks and leave the input sample buffer
* and output array indexes aligned. However, the entire structure (sample buffer, coefficients and output array)
* is inverted (row zero index is M-1 and M-1 is zero index) to facilitate easy loading of data into the input
* sample buffer array using Java's efficient System.arrayCopy() utility.
*
* Note: the filter coefficients array is organized (row,column) for ease of accessing each filter arm row during
* convolution. The samples array is organized (column,row) for ease/efficiency in loading each row of data. Also,
* the input sample array implements a circular buffer where the mCoefficientsMap provides a lookup to each sample
* index based on the current value of the mColumnPointer.
*
* @param coefficients of the prototype filter for this channelizer
*/
private void init(float[] coefficients)
{
int columnCount = (int)Math.ceil((double)coefficients.length / (double)getChannelCount());
int filterRowCount = getChannelCount();
int sampleRowCount = getChannelCount() * 2;
//Sample and filter coefficient arrays are specifically setup as column/row and row/column
mSamples = new float[columnCount][sampleRowCount];
mCoefficients = new float[filterRowCount][columnCount];
mColumnPointer = columnCount - 1;
mBlockSize = getChannelCount();
mRowPointer = mBlockSize;
//Setup the FFT
mFFT = new FloatFFT_1D(getChannelCount());
//Setup the polyphase filters
//Setup the index map for the sample array circular buffer
initIndexMap(columnCount);
int coefficientIndex = 0;
//Load filter columns left-2-right (samples are loaded in reverse column order, right-2-left, for easy convolution)
for(int columnIndex = 0; columnIndex < columnCount; columnIndex++)
{
//Load filter rows bottom to top, opposite of sample row loading and sample output ordering
for(int rowIndex = filterRowCount - 1; rowIndex >= 0; rowIndex--)
{
mCoefficients[rowIndex][columnIndex] = coefficients[coefficientIndex++];
}
}
}
/**
* Generates a circular buffer index map to support lookup of the translated index based on the current column
* pointer and the desired sample index.
*/
private void initIndexMap(int size)
{
mSampleIndexMap = new int[size][size];
for(int row = 0; row < size; row++)
{
for(int column = 0;column < size;column++)
{
int offset = row + column;
mSampleIndexMap[row][column] = ((offset < size) ? offset : offset - size);
}
}
}
public static void main(String[] args)
{
mLog.debug("Starting ...");
int channelCount = 4;
int channelBandwidth = 12500;
int sampleRate = channelCount * channelBandwidth;
int tapsPerFilter = 14;
int channelCenter = (int)(12500 * 1.5);
float[] taps = new float[16];
for(int x = 0; x < taps.length; x++)
{
taps[x] = x;
}
float[] samples = new float[32];
for(int x = 0; x < samples.length; x++)
{
samples[x] = x;
}
try
{
// ComplexPolyphaseChannelizerM2 channelizer = new ComplexPolyphaseChannelizerM2(sampleRate, tapsPerFilter);
ComplexPolyphaseChannelizerM2 channelizer = new ComplexPolyphaseChannelizerM2(taps,sampleRate, channelCount);
// Oscillator oscillator = new Oscillator(channelCenter, 2 * sampleRate);
//
// ComplexBuffer complexBuffer = oscillator.generateComplexBuffer(1000);
ComplexBuffer complexBuffer = new ComplexBuffer(samples);
channelizer.receive(complexBuffer);
}
catch(Exception e)
{
mLog.error("Error", e);
}
mLog.debug("Finished!");
}
}

View File

@@ -0,0 +1,363 @@
/*******************************************************************************
* sdr-trunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
* License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with this program.
* If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.dsp.filter.channelizer.output.IPolyphaseChannelOutputProcessor;
import io.github.dsheirer.dsp.filter.channelizer.output.OneChannelOutputProcessor;
import io.github.dsheirer.dsp.filter.channelizer.output.TwoChannelOutputProcessor;
import io.github.dsheirer.dsp.filter.design.FilterDesignException;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.complex.IComplexBufferProvider;
import io.github.dsheirer.source.Source;
import io.github.dsheirer.source.SourceEvent;
import io.github.dsheirer.source.SourceException;
import io.github.dsheirer.source.tuner.channel.TunerChannel;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.List;
import java.util.concurrent.CopyOnWriteArrayList;
public class PolyphaseChannelManager implements Listener<SourceEvent>
{
private final static Logger mLog = LoggerFactory.getLogger(PolyphaseChannelManager.class);
private static final int CHANNEL_BANDWIDTH = 12500;
private static final double CHANNEL_OVERSAMPLING = 2.0;
private static final int CHANNEL_SAMPLE_RATE = (int)(CHANNEL_BANDWIDTH * CHANNEL_OVERSAMPLING);
private static final int POLYPHASE_FILTER_TAPS_PER_CHANNEL = 17;
private IComplexBufferProvider mComplexBufferProvider;
private List<PolyphaseChannelSource> mChannelSources = new CopyOnWriteArrayList<>();
private ChannelCalculator mChannelCalculator;
private ComplexPolyphaseChannelizerM2 mPolyphaseChannelizer;
private ChannelSourceEventListener mChannelSourceEventListener = new ChannelSourceEventListener();
private ScheduledBufferProcessor mBufferProcessor;
/**
* Polyphase Channel Manager is a DDC channel manager and complex buffer queue/processor for a tuner. This class
* provides DDC polyphase channel sources and wraps a polyphase channelizer processing sample buffers produced by
* the tuner and distributing channelized sample buffers to each allocated DDC polyphase channel source. This
* class is responsible for monitoring the tuner for changes in center frequency and/or sample rate and updating
* active DDC polyphase channel sources accordingly. This class also monitors source event requests and
* notifications received from active DDC polyphase channel sources to adjust sample streams as required.
*
* Note: add this channel manager as a source event listener to the complex buffer provider to ensure this manager
* adapts to changes in source frequency and sample rate.
*
* @param complexBufferProvider that supports register/deregister of this channel manager
* @param frequency of the provided complex buffer samples
* @param sampleRate of the provided complex buffer samples
*/
public PolyphaseChannelManager(IComplexBufferProvider complexBufferProvider, long frequency, int sampleRate)
{
if(complexBufferProvider == null)
{
throw new IllegalArgumentException("Complex buffer provider argument cannot be null");
}
mComplexBufferProvider = complexBufferProvider;
//Ensure that sample rate is a multiple of channel sample rate since polyphase channelizer is M2
if(sampleRate % (CHANNEL_SAMPLE_RATE) != 0)
{
throw new IllegalArgumentException("Tuner sample rate [" + sampleRate + "] must be an even multiple of " +
CHANNEL_SAMPLE_RATE + " Hz");
}
int channelCount = sampleRate / CHANNEL_BANDWIDTH;
mChannelCalculator = new ChannelCalculator(sampleRate, channelCount, frequency, CHANNEL_OVERSAMPLING);
mBufferProcessor = new ScheduledBufferProcessor(500, 100, 50, 50);
updateChannelizer(sampleRate, null);
}
/**
* Provides a Digital Drop Channel (DDC) for the specified tuner channel or returns null if the channel can't be
* sourced due to the current center frequency and/or sample rate.
* @param tunerChannel specifying center frequency and bandwidth.
* @return source or null.
*/
public Source getChannel(TunerChannel tunerChannel)
{
PolyphaseChannelSource channelSource = null;
try
{
List<Integer> polyphaseIndexes = mChannelCalculator.getChannelIndexes(tunerChannel);
IPolyphaseChannelOutputProcessor outputProcessor = getOutputProcessor(polyphaseIndexes);
if(outputProcessor != null)
{
channelSource = new PolyphaseChannelSource(tunerChannel, outputProcessor, mChannelSourceEventListener,
CHANNEL_SAMPLE_RATE);
mChannelSources.add(channelSource);
}
}
catch(IllegalArgumentException iae)
{
mLog.info("Can't provide DDC for " + tunerChannel.toString() + " due to channelizer frequency [" +
mChannelCalculator.getCenterFrequency() + "] and sample rate [" +
(mChannelCalculator.getChannelCount() * mChannelCalculator.getChannelBandwidth()) + "]");
}
return channelSource;
}
/**
* Creates a processor to process the channelizer channel indexes into a composite output stream providing
* channelized complex sample buffers to a registered source listener.
* @param indexes to target by the output processor
* @return output processor compatible with the number of indexes to monitor
*/
private IPolyphaseChannelOutputProcessor getOutputProcessor(List<Integer> indexes)
{
switch(indexes.size())
{
case 1:
return new OneChannelOutputProcessor(CHANNEL_SAMPLE_RATE, indexes);
case 2:
return new TwoChannelOutputProcessor(CHANNEL_SAMPLE_RATE, indexes);
default:
//TODO: create output processor for greater than 2 input channels
mLog.error("Request to create an output processor for unexpected channel index size:" + indexes.size());
return null;
}
}
/**
* Starts/adds the channel source to receive channelized sample buffers, registering with the tuner to receive
* sample buffers when this is the first channel.
*
* @param channelSource to start
*/
private void startChannelSource(PolyphaseChannelSource channelSource)
{
synchronized(mBufferProcessor)
{
mPolyphaseChannelizer.addChannel(channelSource);
if(mPolyphaseChannelizer.getRegisteredChannelCount() == 1)
{
mComplexBufferProvider.addComplexBufferListener(mBufferProcessor);
mBufferProcessor.start();
}
}
}
/**
* Stops/removes the channel source from receiving channelized sample buffers and deregisters from the tuner
* when this is the last channel being sourced.
*
* @param channelSource to stop
*/
private void stopChannelSource(PolyphaseChannelSource channelSource)
{
synchronized(mBufferProcessor)
{
mChannelSources.remove(channelSource);
if(mPolyphaseChannelizer.getRegisteredChannelCount() == 0)
{
mComplexBufferProvider.removeComplexBufferListener(mBufferProcessor);
mBufferProcessor.stop();
}
}
}
/**
* Implements the Listener<SourceEvent> interface for receiving notifications of frequency and sample rate
* changes from the tuner.
*/
@Override
public void receive(SourceEvent sourceEvent)
{
switch(sourceEvent.getEvent())
{
case NOTIFICATION_FREQUENCY_CHANGE:
//Update the index calculator and each of the output processors to update each channel frequency
mChannelCalculator.setCenterFrequency(sourceEvent.getValue().longValue());
updateOutputProcessors(sourceEvent);
break;
case NOTIFICATION_SAMPLE_RATE_CHANGE:
//Update the channelizer (and the output processors) for the new sample rate
updateChannelizer(sourceEvent.getValue().intValue(), sourceEvent);
break;
}
}
/**
* Creates or updates the channelizer to process the incoming sample rate and updates any channel processors.
*
* @param sampleRate to use for the channelizer
* @param sourceEvent containing a notification of sample rate change (optional - can be null). This source
* event will be passed to all output processors to ensure the change is fully promulgated before restarting the
* buffer processor if it was already running.
*/
private void updateChannelizer(int sampleRate, SourceEvent sourceEvent)
{
mChannelCalculator.setSampleRate(sampleRate);
//Lock on the buffer processor to block any channel create/start/stop requests while we configure
synchronized(mBufferProcessor)
{
//Get current processor state so we can restart it after we update the channelizer
boolean bufferProcessorWasRunning = mBufferProcessor.isRunning();
//Stop will be ignored if it is not currently running
mBufferProcessor.stop();
try
{
if(mPolyphaseChannelizer == null)
{
mPolyphaseChannelizer = new ComplexPolyphaseChannelizerM2(sampleRate,
POLYPHASE_FILTER_TAPS_PER_CHANNEL);
mBufferProcessor.setListener(mPolyphaseChannelizer);
}
else
{
mPolyphaseChannelizer.setSampleRate(sampleRate);
}
}
catch(IllegalArgumentException iae)
{
mLog.error("Could not create polyphase channelizer for sample rate [" + sampleRate + "]", iae);
}
catch(FilterDesignException fde)
{
mLog.error("Could not create filter for polyphase channelizer for sample rate [" + sampleRate + "]", fde);
}
if(sourceEvent != null && sourceEvent.getEvent() == SourceEvent.Event.NOTIFICATION_SAMPLE_RATE_CHANGE)
{
updateOutputProcessors(sourceEvent);
}
if(bufferProcessorWasRunning)
{
mBufferProcessor.start();
}
}
}
/**
* Updates each of the output processors for any changes in the tuner's center frequency or sample rate, which
* would cause the output processors to change the polyphase channelizer results channel(s) that the processor is
* consuming
*
* @param sourceEvent that requires an update to the output processors
*/
private void updateOutputProcessors(SourceEvent sourceEvent)
{
for(PolyphaseChannelSource channelSource: mChannelSources)
{
updateOutputProcessor(channelSource);
//Notify the channel source that the frequency or sample rate has changed so that it can reset
//any frequency correction value.
if(sourceEvent != null)
{
try
{
channelSource.process(sourceEvent);
}
catch(SourceException se)
{
mLog.error("Error while notifying polyphase channel source of a source event", se);
}
}
}
}
/**
* Updates the polyphase channel source's output processor due to a change in the center frequency or sample
* rate for the source providing sample buffers to the polyphase channelizer, or whenever the DDC channel's
* center tuned frequency changes.
*
* @param channelSource that requires an update to its output processor
*/
private void updateOutputProcessor(PolyphaseChannelSource channelSource)
{
try
{
List<Integer> indexes = mChannelCalculator.getChannelIndexes(channelSource.getTunerChannel());
long centerFrequency = mChannelCalculator.getCenterFrequencyForIndexes(indexes);
//If the indexes size is the same then update the current processor, otherwise create a new one
if(channelSource.getPolyphaseChannelOutputProcessor().getInputChannelCount() == indexes.size())
{
channelSource.getPolyphaseChannelOutputProcessor().setPolyphaseChannelIndices(indexes);
channelSource.setFrequency(centerFrequency);
}
else
{
channelSource.setPolyphaseChannelOutputProcessor(getOutputProcessor(indexes), centerFrequency);
}
}
catch(IllegalArgumentException iae)
{
mLog.error("Error updating polyphase channel source - can't determine output channel indexes for " +
"updated tuner center frequency and sample rate. Stopping channel source", iae);
stopChannelSource(channelSource);
}
}
/**
* Internal class for handling requests for start/stop sample stream from polyphase channel sources
*/
private class ChannelSourceEventListener implements Listener<SourceEvent>
{
@Override
public void receive(SourceEvent sourceEvent)
{
switch(sourceEvent.getEvent())
{
case REQUEST_START_SAMPLE_STREAM:
if(sourceEvent.hasSource() && sourceEvent.getSource() instanceof PolyphaseChannelSource)
{
startChannelSource((PolyphaseChannelSource)sourceEvent.getSource());
}
else
{
mLog.error("Request to start sample stream for unrecognized source: " +
(sourceEvent.hasSource() ? sourceEvent.getSource().getClass() : "null source"));
}
break;
case REQUEST_STOP_SAMPLE_STREAM:
if(sourceEvent.hasSource() && sourceEvent.getSource() instanceof PolyphaseChannelSource)
{
stopChannelSource((PolyphaseChannelSource)sourceEvent.getSource());
}
else
{
mLog.error("Request to stop sample stream for unrecognized source: " +
(sourceEvent.hasSource() ? sourceEvent.getSource().getClass() : "null source"));
}
break;
default:
mLog.error("Received unrecognized source event from polyphase channel source [" +
sourceEvent.getEvent() + "]");
}
}
}
}

View File

@@ -0,0 +1,222 @@
/*******************************************************************************
* *********************************************************************************************************************
* sdr-trunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
* License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with this program.
* If not, see <http://www.gnu.org/licenses/>
* *********************************************************************************************************************
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.dsp.filter.channelizer.output.IPolyphaseChannelOutputProcessor;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import io.github.dsheirer.sample.complex.ComplexToComplexBufferAssembler;
import io.github.dsheirer.source.SourceEvent;
import io.github.dsheirer.source.tuner.channel.TunerChannel;
import io.github.dsheirer.source.tuner.channel.TunerChannelSource;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class PolyphaseChannelSource extends TunerChannelSource
{
private final static Logger mLog = LoggerFactory.getLogger(PolyphaseChannelSource.class);
private ComplexToComplexBufferAssembler mComplexBufferAssembler = new ComplexToComplexBufferAssembler(2500);
private IPolyphaseChannelOutputProcessor mPolyphaseChannelOutputProcessor;
private double mChannelSampleRate;
private long mCenterFrequency;
private long mFrequencyCorrection;
/**
* Polyphase channelizer tuner channel source implementation. Adapts the channel array output samples from the
* polyphase channelizer into a single channel, or a channel synthesized from two adjacent channels that is
* frequency translated and decimated to a single channel.
*
* @param tunerChannel describing the desired channel frequency and bandwidth/minimum sample rate
* @param outputProcessor - to process polyphase channelizer channel results into a channel stream
* @param producerSourceEventListener to receive source event requests (e.g. start/stop sample stream)
* @param channelSampleRate for the downstream sample output
*/
public PolyphaseChannelSource(TunerChannel tunerChannel, IPolyphaseChannelOutputProcessor outputProcessor,
Listener<SourceEvent> producerSourceEventListener, double channelSampleRate)
{
super(producerSourceEventListener, tunerChannel);
mPolyphaseChannelOutputProcessor = outputProcessor;
mChannelSampleRate = channelSampleRate;
}
/**
* Registers the listener to receive complex sample buffers from this channel source
*/
@Override
public void setListener(final Listener<ComplexBuffer> listener)
{
mComplexBufferAssembler.setListener(listener);
}
/**
* Removes the listener from receiving complex sample buffers from this channel source
*/
@Override
public void removeListener(Listener<ComplexBuffer> listener)
{
mComplexBufferAssembler.setListener(null);
}
/**
* Channel output processor used by this channel source to convert polyphase channel results into a specific
* channel complex buffer output stream.
*/
public IPolyphaseChannelOutputProcessor getPolyphaseChannelOutputProcessor()
{
return mPolyphaseChannelOutputProcessor;
}
/**
* Sets/updates the output processor for this channel source, replacing the existing output processor.
*
* @param outputProcessor to use
* @param frequency center for the channels processed by the output processor
*/
public void setPolyphaseChannelOutputProcessor(IPolyphaseChannelOutputProcessor outputProcessor, long frequency)
{
//Lock on the channel output processor to block the channel results processor thread from servicing the
//channel output processor's buffer while we change out the processors.
synchronized(mPolyphaseChannelOutputProcessor)
{
IPolyphaseChannelOutputProcessor existingProcessor = mPolyphaseChannelOutputProcessor;
//Swap out the processor so that incoming samples can accumulate in the new channel output processor
mPolyphaseChannelOutputProcessor = outputProcessor;
//Remove the buffer overflow listener and clear any previous buffer overflow state
existingProcessor.setSourceOverflowListener(null);
broadcastOverflowState(false);
//Register to receive buffer overflow notifications so we can push them up through the source chain
mPolyphaseChannelOutputProcessor.setSourceOverflowListener(this);
//Fully process the residual channel results buffer of the previous channel output processor
if(existingProcessor != null)
{
existingProcessor.processChannelResults(mComplexBufferAssembler);
}
//Finally, setup the frequency offset for the output processor and reset the frequency correction value
//to allow consumers to adjust and calculate a new frequency correction value, if needed.
mCenterFrequency = frequency;
mFrequencyCorrection = 0;
mPolyphaseChannelOutputProcessor.setFrequencyOffset(getFrequencyOffset());
broadcastConsumerSourceEvent(SourceEvent.frequencyCorrectionChange(mFrequencyCorrection));
}
}
/**
* Primary method for receiving channel results output from a polyphase channelizer. The results array will be
* queued for processing to extract the target channel samples, process them for frequency correction and/or
* channel aggregation, and dispatch the results to the downstream sample listener/consumer.
*
* @param channelResultsBuffer containing the polyphase channelizer output for a single complex sample period.
*/
public void receiveChannelResults(float[] channelResultsBuffer)
{
mPolyphaseChannelOutputProcessor.receiveChannelResults(channelResultsBuffer);
}
/**
* Downstream channel sample rate
*
* @return sample rate in Hertz
*/
@Override
public double getSampleRate()
{
return mChannelSampleRate;
}
/**
* Sets the downstream channel sample rate
* @param sampleRate in hertz
*/
@Override
protected void setSampleRate(double sampleRate)
{
mChannelSampleRate = sampleRate;
}
@Override
public void dispose()
{
mLog.info("Polyphase channel source received request to dispose");
}
@Override
public long getFrequencyCorrection()
{
return mFrequencyCorrection;
}
/**
* Adjusts the frequency correction value that is being applied to the channelized output stream by the
* polyphase channel output processor.
*
* @param value to apply for frequency correction in hertz
*/
protected void setFrequencyCorrection(long value)
{
mFrequencyCorrection = value;
updateFrequencyOffset();
broadcastConsumerSourceEvent(SourceEvent.frequencyCorrectionChange(mFrequencyCorrection));
}
/**
* Sets the center frequency for the incoming sample stream channel results and resets frequency correction to zero.
* @param frequency in hertz
*/
@Override
protected void setFrequency(long frequency)
{
mCenterFrequency = frequency;
setFrequencyCorrection(0);
}
@Override
protected void processSamples()
{
//Lock on the output processor so that it can't be changed out in the middle of processing
synchronized(mPolyphaseChannelOutputProcessor)
{
mPolyphaseChannelOutputProcessor.processChannelResults(mComplexBufferAssembler);
}
}
/**
* Calculates the frequency offset required to mix the incoming signal to the desired output baseband
*/
private long getFrequencyOffset()
{
return getTunerChannel().getFrequency() - mCenterFrequency + mFrequencyCorrection;
}
/**
* Updates the frequency offset being applied by the output processor. Calculates the offset using the center
* frequency of the incoming polyphase channel results and the desired output channel frequency adjusted by any
* requested frequency correction value.
*/
private void updateFrequencyOffset()
{
synchronized(mPolyphaseChannelOutputProcessor)
{
mPolyphaseChannelOutputProcessor.setFrequencyOffset(getFrequencyOffset());
}
}
}

View File

@@ -0,0 +1,168 @@
/*******************************************************************************
* sdr-trunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
* License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along with this program.
* If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.OverflowableTransferQueue;
import io.github.dsheirer.sample.real.IOverflowListener;
import io.github.dsheirer.util.ThreadPool;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.ScheduledFuture;
import java.util.concurrent.TimeUnit;
import java.util.concurrent.atomic.AtomicBoolean;
public class ScheduledBufferProcessor<E> implements Listener<E>
{
private final static Logger mLog = LoggerFactory.getLogger(ScheduledBufferProcessor.class);
private OverflowableTransferQueue<E> mQueue;
private Listener<E> mListener;
private ScheduledFuture<?> mScheduledFuture;
private long mDistributionInterval;
private int mMaxBuffersPerInterval;
private AtomicBoolean mRunning = new AtomicBoolean();
/**
* Scheduled Buffer Processor combines an internal overflowable buffer with a scheduled runnable processing task
* for periodically distributing internally queued elements to the registered listener. This processor provides
* a convenient way to create a thread-safe buffer for receiving elements from one thread/runnable and then
* distributing those elements to a registered listener where distribution occurs on a separate scheduled thread
* pool runnable thread. This allows the calling input thread to quickly return without incurring any subsequent
* processing workload.
*
* The internal queue is an overflowable queue implementation that allows a listener to be registered to receive
* notifications of overflow and reset state. Queue sizing parameters are specified in the constructor.
*
* @param maximumSize of the internal queue (overflow happens when this is exceeded)
* @param resetThreshold of the internal queue (overflow reset happens once queue size falls below this threshold
* @param distributionInterval in milliseconds specifying how often the internal scheduled thread pool buffer
* distributor should run
* @param maxBuffersPerInterval maximum number of buffers to pull from the internal queue during each distribution
* interval
*/
public ScheduledBufferProcessor(int maximumSize, int resetThreshold, long distributionInterval, int maxBuffersPerInterval)
{
mQueue = new OverflowableTransferQueue<>(maximumSize, resetThreshold);
mQueue.setOverflowListener(new IOverflowListener()
{
@Override
public void sourceOverflow(boolean overflow)
{
mLog.error("Overflow state changed - overflow:" + overflow);
}
});
mDistributionInterval = distributionInterval;
mMaxBuffersPerInterval = maxBuffersPerInterval;
}
/**
* Sets the listener to receive notifications of buffer overflow and/or reset. Note: this method can also be
* invoked with a null argument to clear the previously registered overflow listener.
*
* @param listener to receive overflow/reset notifications.
*/
public void setOverflowListener(IOverflowListener listener)
{
mQueue.setOverflowListener(listener);
}
/**
* Sets or changes the listener to receive buffers from this processor.
* @param listener to receive buffers
*/
public void setListener(Listener<E> listener)
{
mListener = listener;
}
/**
* Primary input method for adding buffers to this processor. Note: incoming buffers will be ignored if this
* processor is in a stopped state. You must invoke start() to allow incoming buffers and initiate buffer
* processing.
*
* @param e to enqueue for distribution to a registered listener
*/
@Override
public void receive(E e)
{
mQueue.offer(e);
}
/**
* Starts this buffer processor and allows queuing of incoming buffers.
*/
public void start()
{
if(mRunning.compareAndSet(false, true))
{
mScheduledFuture = ThreadPool.SCHEDULED.scheduleAtFixedRate(new Processor(), 0, mDistributionInterval, TimeUnit.MILLISECONDS);
}
}
/**
* Stops this buffer processor, stops queuing of incoming buffers and clears the internal buffer queue.
*/
public void stop()
{
if(mRunning.compareAndSet(true, false))
{
if(mScheduledFuture != null)
{
mScheduledFuture.cancel(true);
mScheduledFuture = null;
}
mQueue.clear();
}
}
/**
* Indicates if this processor is currently running
*/
public boolean isRunning()
{
return mRunning.get();
}
/**
* Processor to service the buffer queue and distribute the buffers to the registered listener
*/
class Processor implements Runnable
{
private List<E> mBuffers = new ArrayList<>();
@Override
public void run()
{
mQueue.drainTo(mBuffers, mMaxBuffersPerInterval);
if(mListener != null)
{
for(E buffer: mBuffers)
{
mListener.receive(buffer);
}
}
mBuffers.clear();
}
}
}

View File

@@ -0,0 +1,170 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer.output;
import io.github.dsheirer.dsp.mixer.Oscillator;
import io.github.dsheirer.sample.OverflowableTransferQueue;
import io.github.dsheirer.sample.complex.Complex;
import io.github.dsheirer.sample.complex.ComplexSampleListener;
import io.github.dsheirer.sample.real.IOverflowListener;
import io.github.dsheirer.source.Source;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.ArrayList;
import java.util.List;
public abstract class ChannelOutputProcessor implements IPolyphaseChannelOutputProcessor
{
private final static Logger mLog = LoggerFactory.getLogger(ChannelOutputProcessor.class);
private OverflowableTransferQueue<float[]> mChannelResultsQueue;
private List<float[]> mChannelResultsToProcess = new ArrayList<>();
private int mMaxResultsToProcess;
private int mInputChannelCount;
//TODO: swap this out and use the LowPhaseNoiseOscillator
private Oscillator mFrequencyCorrectionMixer;
private boolean mFrequencyCorrectionEnabled;
/**
* Base class for polyphase channelizer output channel processing. Provides built-in frequency translation
* oscillator support to apply frequency correction to the channel sample stream as requested by sample consumer.
*
* @param inputChannelCount is the number of input channels for this output processor
* @param sampleRate of the output channel. This is used to match the oscillator's sample rate to the output
* channel sample rate for frequency translation/correction.
*/
public ChannelOutputProcessor(int inputChannelCount, double sampleRate)
{
mInputChannelCount = inputChannelCount;
mFrequencyCorrectionMixer = new Oscillator(0, sampleRate);
mMaxResultsToProcess = (int)(sampleRate / 10) * 2; //process at 100 millis interval, twice the expected inflow rate
mChannelResultsQueue = new OverflowableTransferQueue<>((int)(sampleRate * 3), (int)(sampleRate * 0.5));
}
/**
* Sets the frequency offset to apply to the incoming samples to mix the desired signal to baseband.
* @param frequencyOffset in hertz
*/
@Override
public void setFrequencyOffset(long frequencyOffset)
{
mFrequencyCorrectionMixer.setFrequency(frequencyOffset);
mFrequencyCorrectionEnabled = (frequencyOffset != 0);
}
@Override
public void receiveChannelResults(float[] channelResults)
{
mChannelResultsQueue.offer(channelResults);
}
/**
* Processes all enqueued polyphase channelizer results until the internal queue is empty
* @param listener to receive the processed channel results
*/
@Override
public void processChannelResults(ComplexSampleListener listener)
{
try
{
int count = mChannelResultsQueue.drainTo(mChannelResultsToProcess, mMaxResultsToProcess);
if(count > 0)
{
process(mChannelResultsToProcess, listener);
}
}
catch(Throwable throwable)
{
mLog.error("Error while processing polyphase channel samples", throwable);
}
mChannelResultsToProcess.clear();
}
/**
* Sub-class implementation to process one polyphase channelizer result array.
* @param channelResults to process
* @param listener to receive the channel complex sample
*/
public abstract void process(List<float[]> channelResults, ComplexSampleListener listener);
/**
* Sets the overflow listener to monitor the internal channelizer channel results queue overflow state
*/
public void setSourceOverflowListener(Source source)
{
mChannelResultsQueue.setSourceOverflowListener(source);
}
/**
* Removes the overflow listener from monitoring the internal channelizer channel results queue overflow state
*/
public void removeOverflowListener(IOverflowListener listener)
{
mChannelResultsQueue.setOverflowListener(null);
}
@Override
public int getInputChannelCount()
{
return mInputChannelCount;
}
/**
* Provides frequency correction of the inphase component or simply returns the inphase argument value when no
* frequency correction value has been specified.
*
* @param inphase value of the sample to translate/correct
* @param quadrature value of the sample to translate/correct
* @return corrected inphase value
*/
protected float getFrequencyCorrectedInphase(float inphase, float quadrature)
{
return Complex.multiplyInphase(inphase, quadrature, mFrequencyCorrectionMixer.inphase(),
mFrequencyCorrectionMixer.quadrature());
}
/**
* Provides frequency correction of the quadrature component or simply returns the quadrature argument value when no
* frequency correction value has been applied,
*
* @param inphase value of the sample to translate/correct
* @param quadrature value of the sample to translate/correct
* @return translated inphase value
*/
protected float getFrequencyCorrectedQuadrature(float inphase, float quadrature)
{
return Complex.multiplyQuadrature(inphase, quadrature, mFrequencyCorrectionMixer.inphase(),
mFrequencyCorrectionMixer.quadrature());
}
/**
* Indicates if this channel has a frequency correction offset
*/
public boolean hasFrequencyCorrection()
{
return mFrequencyCorrectionEnabled;
}
}

View File

@@ -0,0 +1,66 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer.output;
import io.github.dsheirer.sample.complex.ComplexSampleListener;
import io.github.dsheirer.source.Source;
import java.util.List;
public interface IPolyphaseChannelOutputProcessor
{
/**
* Receive and enqueue output results from the polyphase analysis channelizer
* @param channelResults to enqueue
*/
void receiveChannelResults(float[] channelResults);
/**
* Process the channel output channel results queue and deliver the output to the listener
*
* @param listener to receive the processed channel results
*/
void processChannelResults(ComplexSampleListener listener);
/**
* Sets the desired frequency offset from center. The samples will be mixed with an oscillator set to this offset
* frequency to produce an output where the desired signal is centered in the passband.
*
* @param frequency in hertz
*/
void setFrequencyOffset(long frequency);
/**
* Indicates the number of input channels processed by this output processor
* @return
*/
int getInputChannelCount();
/**
* Updates the input polyphase channel index(es) used by this output processor
* @param indexes
*/
void setPolyphaseChannelIndices(List<Integer> indexes);
/**
* Sets the listener to receive buffer overflow notifications
* @param source to receive overflow notifications
*/
void setSourceOverflowListener(Source source);
}

View File

@@ -0,0 +1,96 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer.output;
import io.github.dsheirer.sample.complex.ComplexSampleListener;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.List;
public class OneChannelOutputProcessor extends ChannelOutputProcessor
{
private final static Logger mLog = LoggerFactory.getLogger(OneChannelOutputProcessor.class);
private int mChannelOffset;
/**
* Processor to extract a single channel from a polyphase channelizer and produce an output I/Q complex sample
* from each polyphase channelizer output results array.
*
* @param sampleRate of the output sample stream.
* @param channelIndexes containing a single channel index.
*/
public OneChannelOutputProcessor(double sampleRate, List<Integer> channelIndexes)
{
super(1, sampleRate);
setPolyphaseChannelIndices(channelIndexes);
mLog.debug("Sample Rate:" + sampleRate + " PolyPhase Channel Indexes:" + channelIndexes);
}
/**
* Updates this processor to extract a single, specified channel index.
* @param indexes containing a single channel index value.
* @throws IllegalArgumentException if the list of indexes does not contain a single channel index.
*/
public void setPolyphaseChannelIndices(List<Integer> indexes)
{
if(indexes.size() != 1)
{
throw new IllegalArgumentException("Single channel output processor requires a single index to " +
"process - provided indexes " + indexes.toString());
}
//Set the channelized output results offset to twice the channel index to account for each channel having
//an I/Q pair
mChannelOffset = indexes.get(0) * 2;
}
/**
* Extract the channel from the channel results array, apply frequency translation, and deliver the
* extracted frequency-corrected channel I/Q sample set to the complex sample listener.
*
* @param channelResults to process containing a list of channel array of I/Q sample pairs (I0,Q0,I1,Q1...In,Qn)
* @param listener to receive the extracted, frequency-translated channel results
*/
@Override
public void process(List<float[]> channelResults, ComplexSampleListener listener)
{
for(float[] channels: channelResults)
{
if(channels.length < mChannelOffset + 1)
{
throw new IllegalArgumentException("Polyphase channelizer output channels array is not large enough to " +
"cover this single channel output processor with channel offset [" + mChannelOffset + "]");
}
float i = channels[mChannelOffset];
float q = channels[mChannelOffset + 1];
if(hasFrequencyCorrection())
{
listener.receive(getFrequencyCorrectedInphase(i, q), getFrequencyCorrectedQuadrature(i, q));
}
else
{
listener.receive(i, q);
}
}
}
}

View File

@@ -0,0 +1,91 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.channelizer.output;
import io.github.dsheirer.sample.complex.ComplexSampleListener;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.List;
public class TwoChannelOutputProcessor extends ChannelOutputProcessor
{
private final static Logger mLog = LoggerFactory.getLogger(TwoChannelOutputProcessor.class);
private int mChannelOffset1;
private int mChannelOffset2;
/**
* Processor to extract two channels from a polyphase channelizer, synthesize/recombine the channels, apply
* frequency translation and frequency correction, down-sample by a factor of two and output an I/Q complex sample.
*
* @param sampleRate of the output sample stream.
* @param channelIndexes containing two channel indices.
*/
public TwoChannelOutputProcessor(int sampleRate, List<Integer> channelIndexes)
{
//Set the frequency correction oscillator to 2 x output sample rate since we'll be correcting the frequency
//after synthesizing both input channels
super(2, sampleRate * 2);
setPolyphaseChannelIndices(channelIndexes);
}
/**
* Updates this processor to extract the two specified channel indexes.
*
* @param indexes containing a single channel index value.
* @throws IllegalArgumentException if the list of indexes does not contain a single channel index.
*/
public void setPolyphaseChannelIndices(List<Integer> indexes)
{
if(indexes.size() != 2)
{
throw new IllegalArgumentException("Double channel output processor requires two indexes to " +
"process - provided indexes " + indexes.toString());
}
//Set the channelized output results offsets to twice the channel index to account for each channel having
//an I/Q pair
mChannelOffset1 = indexes.get(0) * 2;
mChannelOffset2 = indexes.get(1) * 2;
}
/**
* Extract the channel from the channel results array, apply frequency translation, and deliver the
* extracted frequency-corrected channel I/Q sample set to the complex sample listener.
*
* @param channelResults to process containing an array of channel I/Q sample pairs (I0,Q0,I1,Q1...In,Qn)
* @param listener to receive the extracted, frequency-translated channel results
*/
@Override
public void process(List<float[]> channelResults, ComplexSampleListener listener)
{
mLog.debug("I got channel results - doing nothing !!");
// if(channels.length < mChannelOffset + 1)
// {
// throw new IllegalArgumentException("Polyphase channelizer output channels array is not large enough to " +
// "cover this single channel output processor with channel offset [" + mChannelOffset + "]");
// }
//
// float i = channels[mChannelOffset];
// float q = channels[mChannelOffset + 1];
// listener.receive(getFrequencyCorrectedInphase(i, q), getFrequencyCorrectedQuadrature(i, q));
}
}

View File

@@ -0,0 +1,51 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.correction;
import io.github.dsheirer.sample.complex.ComplexBuffer;
public class IQCorrectionFilter
{
private float mRatio = 0.00001f;
private float mAverageInphase = 0.0f;
private float mAverageQuadrature = 0.0f;
public IQCorrectionFilter(float ratio)
{
mRatio = ratio;
}
public IQCorrectionFilter()
{
}
public void filter(ComplexBuffer buffer)
{
float[] samples = buffer.getSamples();
for(int x = 0; x < samples.length; x += 2)
{
mAverageInphase = mAverageInphase + (mRatio * (samples[x] - mAverageInphase));
samples[x] -= mAverageInphase;
mAverageQuadrature = mAverageQuadrature + (mRatio * (samples[x + 1] - mAverageQuadrature));
samples[x + 1] -= mAverageQuadrature;
}
}
}

View File

@@ -0,0 +1,37 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.fir.remez;
import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification;
public class PolyphaseChannelizerDesigner extends FIRFilterSpecification
{
/**
* Constructs a Polyphase TunerChannelizer Remez filter specification.
*
* Constructor is package private -- use the PolyphaseChannelizerFilterFactory to create a filter
*
* @param order of the filter
* @param gridDensity to use when calculating the coefficients
*/
PolyphaseChannelizerDesigner(int order, int gridDensity)
{
super(FIRLinearPhaseFilterType.TYPE_1_ODD_LENGTH_EVEN_ORDER_SYMMETRICAL, order, gridDensity);
}
}

View File

@@ -0,0 +1,97 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.fir.remez;
import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class PolyphaseChannelizerFilterFactory
{
private final static Logger mLog = LoggerFactory.getLogger(PolyphaseChannelizerFilterFactory.class);
private static final double OBJECTIVE_BAND_EDGE_COEFFICIENT_AMPLITUDE = Math.sqrt(2.0) / 2.0; //.707xxx
public static float[] getFilter(int sampleRate, int channelBandwidth, double alpha)
{
FIRLinearPhaseFilterType type = FIRLinearPhaseFilterType.TYPE_2_EVEN_LENGTH_ODD_ORDER_SYMMETRICAL;
int passBandStart = 0;
int gridDensity = 60;
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 = channelBandwidth + rolloffFrequency;
double passRipple = 0.01d;
double transitionRipple = 0.02d;
double stopRipple = 0.001d; //Approximately 90 dB of attenuation
int order = FIRFilterSpecification.estimateFilterOrder( sampleRate, passBandStop, stopBandStart,
passRipple, stopRipple );
// order = 1001;
// Ensure odd order since we're designing a Type 2 filter
if(order % 2 == 0)
{
order++;
}
mLog.info("Filter Order: " + order);
FIRFilterSpecification specification = new PolyphaseChannelizerDesigner(order, gridDensity);
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, transitionRipple, 5.0);
FIRFilterSpecification.FrequencyBand stopBand = new FIRFilterSpecification.FrequencyBand(sampleRate,
stopBandStart, (int)(sampleRate / 2), 0.0, stopRipple);
specification.addFrequencyBand(passBand);
specification.addFrequencyBand(transitionBand);
specification.addFrequencyBand(stopBand);
double bandEdgeFrequency = Math.cos(Math.PI * (double)channelBandwidth / (double)(sampleRate / 2));
float[] filter = null;
try
{
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);
}
catch(Exception e)
{
mLog.error("Error designing filter", e);
}
return filter;
}
}

View File

@@ -211,15 +211,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
{

View File

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

View File

@@ -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 <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.filter.fir.remez;
import io.github.dsheirer.dsp.filter.design.FilterDesignException;
import io.github.dsheirer.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<Integer> 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<Integer> getInitialExtremalIndices()
{
int count = mSpecification.getExtremaCount();
int density = mSpecification.getGridDensity();
List<Integer> extremalIndices = new ArrayList<>();
for(int i = 0; i < count; i++)
{
extremalIndices.add(i * density);
}
return extremalIndices;
}
/**
* Calculates the actual frequency response of a dense frequency grid using the polynomial
* represented by the extremal indices.
*/
private void calculateGridFrequencyResponse()
{
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<Integer> indicesToRemove = new ArrayList<>();
Iterator<Integer> it = mExtremalIndices.iterator();
Integer current = it.next();
Integer next;
boolean positiveAxis = mGridErrors[current] > 0.0;
while(it.hasNext())
{
next = it.next();
//Check for consecutive (redundant) indices on same side of zero axis - retain largest
if(!(positiveAxis ^ (mGridErrors[next] > 0.0)))
{
if(Math.abs(mGridErrors[next]) <= Math.abs(mGridErrors[current]))
{
//Remove next error index that is less than current index and on
it.remove();
next = current;
}
else
{
//Since we can't remove the current index with the iterator, queue for later removal
indicesToRemove.add(current);
}
}
else
{
positiveAxis = !positiveAxis;
}
current = next;
}
//Remove redundant extremal indices that couldn't be removed via the iterator
mExtremalIndices.removeAll(indicesToRemove);
//Truncate the list to one larger than needed by removing excess tailing indices
while(mExtremalIndices.size() > mSpecification.getExtremaCount())
{
mExtremalIndices.remove(mExtremalIndices.size() - 1);
}
//If we have one too many indices, delete the smaller of the first or last
if(mExtremalIndices.size() > mSpecification.getExtremaCount())
{
int lastIndex = mExtremalIndices.size() - 1;
if(Math.abs(mGridErrors[mExtremalIndices.get(0)]) >
Math.abs(mGridErrors[mExtremalIndices.get(lastIndex)]))
{
mExtremalIndices.remove(lastIndex);
}
else
{
mExtremalIndices.remove(0);
}
}
//Detect if we have too few indices
if(mExtremalIndices.size() < mSpecification.getExtremaCount())
{
throw new FilterDesignException("Couldn't find the minimum extremal frequencies in "
+ "error set - 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;
}
}

View File

@@ -142,7 +142,7 @@ public class AutomaticGainControl
public float process( float currentSample )
{
float delayedSample = mDelayBuffer.get( currentSample );
float delayedSample = mDelayBuffer.putAndGet( currentSample );
double gain = MANUAL_AGC_GAIN;

View File

@@ -74,7 +74,7 @@ public class ComplexFeedForwardGainControl implements ComplexSampleListener
}
/* Replace oldest envelope value with current envelope value */
float oldestEnvelope = mEnvelopeHistory.get( envelope );
float oldestEnvelope = mEnvelopeHistory.putAndGet( envelope );
/* If the oldest envelope value was the max envelope value, then we
* have to rediscover the max value from the envelope history */

View File

@@ -0,0 +1,188 @@
/*******************************************************************************
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.mixer;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class LowPhaseNoiseOscillator
{
private final static Logger mLog = LoggerFactory.getLogger(LowPhaseNoiseOscillator.class);
private static final double TWO_PI = 2.0 * Math.PI;
private static final float THREE_HALVES = 3.0f / 2.0f;
private int mSampleRate;
private long mFrequency;
private float mInphase = 1.0f;
private float mQuadrature = 0.0f;
private float mPreviousInphase = 1.0f;
private float mPreviousQuadrature = 0.0f;
private float mCosineAngle;
private float mSineAngle;
private float mGain = 1.0f;
/**
* Ultra-low phase noise complex oscillator as described in Digital Signal Processing 3e, Lyons, p.786
*
* @param sampleRate for the oscillator
* @param frequency to generate
*/
public LowPhaseNoiseOscillator(int sampleRate, long frequency)
{
mSampleRate = sampleRate;
mFrequency = frequency;
update();
}
/**
* Sets the sample rate for this oscillator
*/
public void setSampleRate(int sampleRate)
{
mSampleRate = sampleRate;
update();
}
/**
* Sample rate for this oscillator
*/
public int getSampleRate()
{
return mSampleRate;
}
/**
* Sets the tone frequency generated by this oscillator
*/
public void setFrequency(long frequency)
{
mFrequency = frequency;
update();
}
/**
* Current tone frequency for this oscillator
*/
public long getFrequency()
{
return mFrequency;
}
/**
* Updates the internal values after a frequency or sample rate change
*/
private void update()
{
double anglePerSample = TWO_PI * mFrequency / mSampleRate;
mCosineAngle = (float)Math.cos(anglePerSample);
mSineAngle = (float)Math.sin(anglePerSample);
}
private void rotate()
{
mInphase = ((mPreviousInphase * mCosineAngle) - (mPreviousQuadrature * mSineAngle)) * mGain;
mQuadrature = ((mPreviousInphase * mSineAngle) + (mPreviousQuadrature * mCosineAngle)) * mGain;
mPreviousInphase = mInphase;
mPreviousQuadrature = mQuadrature;
mGain = THREE_HALVES - ((mPreviousInphase * mPreviousInphase) + (mPreviousQuadrature * mPreviousQuadrature));
}
/**
* Generates an array of samples with length of 2 x sample count containing interleaved inphase and quadrature
* samples.
*
* @param sampleCount number of samples to generate
* @return array that is sampleCount * 2 in length
*/
public float[] generate(int sampleCount)
{
float[] samples = new float[sampleCount * 2];
for(int x = 0; x < sampleCount; x++)
{
int index = x * 2;
samples[index] = mInphase;
samples[index + 1] = mQuadrature;
rotate();
}
return samples;
}
public static void main(String[] args)
{
int sampleRate = 10000000;
int frequency = 24000;
int sampleSize = 65536;
int warmUpIterations = 100000;
int testIterations = 500000;
LowPhaseNoiseOscillator o = new LowPhaseNoiseOscillator(sampleRate, frequency);
mLog.debug("Warming up ...");
for(int x = 0; x < warmUpIterations; x++)
{
float[] samples = o.generate(sampleSize);
}
mLog.debug("Testing ...");
long start = System.nanoTime();
for(int x = 0; x < testIterations; x++)
{
float[] samples = o.generate(sampleSize);
}
long duration = System.nanoTime() - start;
double durationSeconds = (double)duration / 1000000000.0;
mLog.debug("New Took:" + durationSeconds + " seconds to generate " + (testIterations * sampleSize) + " samples");
Oscillator o2 = new Oscillator(frequency, sampleRate);
mLog.debug("Warming up ...");
for(int x = 0; x < warmUpIterations; x++)
{
ComplexBuffer samples = o2.generateComplexBuffer(sampleSize);
}
mLog.debug("Testing ...");
long start2 = System.nanoTime();
for(int x = 0; x < testIterations; x++)
{
float[] samples = o.generate(sampleSize);
}
long duration2 = System.nanoTime() - start2;
double durationSeconds2 = (double)duration2 / 1000000000.0;
mLog.debug("New Took:" + durationSeconds2 + " seconds to generate " + (testIterations * sampleSize) + " samples");
}
}

View File

@@ -1,112 +1,152 @@
/*******************************************************************************
* SDR Trunk
* Copyright (C) 2014 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
* sdrtrunk
* Copyright (C) 2014-2017 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*
******************************************************************************/
package io.github.dsheirer.dsp.mixer;
import io.github.dsheirer.sample.complex.Complex;
import io.github.dsheirer.sample.complex.ComplexBuffer;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
public class Oscillator
{
private final static Logger mLog = LoggerFactory.getLogger( Oscillator.class );
private final static Logger mLog = LoggerFactory.getLogger(Oscillator.class);
private double mFrequency;
private double mSampleRate;
private boolean mEnabled = true;
private double mFrequency;
private double mSampleRate;
private Complex mAnglePerSample;
private Complex mCurrentAngle = new Complex( 0.0f, -1.0f );
private Complex mAnglePerSample;
private Complex mCurrentAngle = new Complex(0.0f, -1.0f);
/**
* Oscillator produces complex or float samples corresponding to a sine wave
* oscillating at the specified frequency and sample rate
*
* @param frequency - positive or negative frequency in hertz
* @param sampleRate - in hertz
*/
public Oscillator( long frequency, int sampleRate )
{
mSampleRate = (double)sampleRate;
mFrequency = (double)frequency;
update();
}
/**
* Updates the internal values after a frequency or sample rate change
*/
private void update()
{
float anglePerSample =
(float)( 2.0d * Math.PI * mFrequency / mSampleRate );
/**
* Oscillator produces complex or float samples corresponding to a sine wave
* oscillating at the specified frequency and sample rate
*
* @param frequency - positive or negative frequency in hertz
* @param sampleRate - in hertz
*/
public Oscillator(long frequency, double sampleRate)
{
mEnabled = frequency != 0;
mAnglePerSample = Complex.fromAngle( anglePerSample );
}
mSampleRate = sampleRate;
mFrequency = (double)frequency;
/**
* Sets or changes the frequency of this oscillator
*/
public void setFrequency( long frequency )
{
mFrequency = (double)frequency;
update();
}
update();
}
/**
* Sets or changes the sample rate of this oscillator
*/
public void setSampleRate( int sampleRate )
{
mSampleRate = (double)sampleRate;
update();
}
/**
* Indicates if this oscillator has a non-zero frequency setting. When used as a frequency translating
* mixer, setting a frequency of zero indicates that no frequency correction is needed, so this flag can
* indicate when to apply the oscillator to frequency correct samples.
*
* @return true if the current frequency value is non-zero.
*/
public boolean enabled()
{
return mEnabled;
}
/**
* Steps the current angle by the angle per sample amount
*/
public void rotate()
{
mCurrentAngle.multiply( mAnglePerSample );
mCurrentAngle.fastNormalize();
}
/**
* Updates the internal values after a frequency or sample rate change
*/
private void update()
{
float anglePerSample =
(float)(2.0d * Math.PI * mFrequency / mSampleRate);
public float inphase()
{
return mCurrentAngle.inphase();
}
public float quadrature()
{
return mCurrentAngle.quadrature();
}
/**
* Get next complex sample
*/
public Complex getComplex()
{
return mCurrentAngle.copy();
}
mAnglePerSample = Complex.fromAngle(anglePerSample);
}
/**
* Get the next float sample
*/
public float getFloat()
{
return mCurrentAngle.real();
}
/**
* Sets or changes the frequency of this oscillator
*/
public void setFrequency(long frequency)
{
mEnabled = frequency != 0;
mFrequency = (double)frequency;
update();
}
/**
* Sets or changes the sample rate of this oscillator
*/
public void setSampleRate(double sampleRate)
{
mSampleRate = sampleRate;
update();
}
/**
* Steps the current angle by the angle per sample amount
*/
public void rotate()
{
mCurrentAngle.multiply(mAnglePerSample);
mCurrentAngle.fastNormalize();
}
public float inphase()
{
return mCurrentAngle.inphase();
}
public float quadrature()
{
return mCurrentAngle.quadrature();
}
/**
* Get next complex sample
*/
public Complex getComplex()
{
return mCurrentAngle.copy();
}
/**
* Get the next float sample
*/
public float getFloat()
{
return mCurrentAngle.real();
}
/**
* Generates the requested number of complex samples and returns them in a complex buffer
*
* @param sampleCount to generate
* @return complex buffer containing the requested complex samples
*/
public ComplexBuffer generateComplexBuffer(int sampleCount)
{
float[] buffer = new float[sampleCount * 2];
for(int x = 0; x < sampleCount; x++)
{
int index = 2 * x;
buffer[index] = inphase();
buffer[index + 1] = quadrature();
rotate();
}
return new ComplexBuffer(buffer);
}
}

View File

@@ -41,7 +41,7 @@ import io.github.dsheirer.instrument.tap.stream.EyeDiagramDataTap;
import io.github.dsheirer.sample.Listener;
import io.github.dsheirer.sample.complex.Complex;
import io.github.dsheirer.sample.complex.ComplexSampleListener;
import io.github.dsheirer.source.tuner.frequency.IFrequencyChangeListener;
import io.github.dsheirer.source.ISourceEventListener;
import java.util.ArrayList;
import java.util.Arrays;
@@ -72,7 +72,7 @@ public class LSMDemodulator implements Instrumentable, ComplexSampleListener
private CostasLoop mCostasLoop = new CostasLoop();
private IFrequencyChangeListener mFrequencyChangeListener;
private ISourceEventListener mFrequencyChangeListener;
private EyeDiagramDataTap mEyeDiagramDataTap;
@@ -146,12 +146,12 @@ public class LSMDemodulator implements Instrumentable, ComplexSampleListener
}
}
public void addListener( IFrequencyChangeListener listener )
public void addListener( ISourceEventListener listener )
{
mFrequencyChangeListener = listener;
}
public void removeListener( IFrequencyChangeListener listener )
public void removeListener( ISourceEventListener listener )
{
mFrequencyChangeListener = null;
}
@@ -201,8 +201,8 @@ public class LSMDemodulator implements Instrumentable, ComplexSampleListener
// if( correctionNeeded && mFrequencyChangeListener != null )
// {
// mLog.debug( "Issuing Frequency Correction: " + mFrequencyError );
// mFrequencyChangeListener.frequencyChanged(
// new FrequencyChangeEvent( Attribute.FREQUENCY_ERROR,
// mFrequencyChangeListener.process(
// new SourceEvent( Attribute.FREQUENCY_ERROR,
// mFrequencyError ) );
// }
//