Initial code check in continued

This commit is contained in:
dennis.sheirer
2014-04-26 19:48:24 +00:00
parent 960a6f9952
commit 3dfedcf8d1
474 changed files with 76613 additions and 0 deletions

View File

@@ -0,0 +1,117 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
import sample.Listener;
public class BooleanZeroCrossingDetector implements Listener<Boolean>
{
private int mBaudLength;
private int mCounter = 0;
private int mThreshold;
private Listener<Boolean> mListener;
private Cycle mCycle = Cycle.POSITIVE;
/**
* Implements a zero crossing detector for use with squared, filtered
* binary values. Makes a baud decision at 100% of the baud
* length, or at 80% of the baud length during a zero crossing. Counter
* is reset at each zero crossing.
*
* @param baudLength - estimated number of boolean samples per baud
*/
public BooleanZeroCrossingDetector( int baudLength )
{
mBaudLength = baudLength;
mThreshold = (int)( baudLength * .80 );
}
@Override
public void receive( Boolean sample )
{
mCounter++;
switch( mCycle )
{
case POSITIVE:
if( sample )
{
if( mCounter >= mBaudLength )
{
send( true );
mCounter -= mBaudLength;
}
}
else
{
if( mCounter >= mThreshold )
{
send( true );
mCounter -= mBaudLength;
}
mCycle = Cycle.NEGATIVE;
}
break;
case NEGATIVE:
if( sample )
{
if( mCounter >= mThreshold )
{
send( false );
mCounter -= mBaudLength;
}
mCycle = Cycle.POSITIVE;
}
else
{
if( mCounter >= mBaudLength )
{
send( false );
mCounter -= mBaudLength;
}
}
break;
}
}
private void send( boolean bit )
{
if( mListener != null )
{
mListener.receive( bit );
}
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void removeListener( Listener<Boolean> listener )
{
mListener = null;
}
public enum Cycle { POSITIVE, NEGATIVE };
}

View File

@@ -0,0 +1,69 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
import sample.Listener;
import buffer.FloatCircularAveragingBuffer;
public class FloatToBinarySlicer implements Listener<Float>
{
private FloatCircularAveragingBuffer mAveragingBuffer;
private Listener<Boolean> mListener;
/**
* Provides binary slicing of the stream of floats, using a circular buffer
* to calculate a running average, then basing the slice decision based
* on the oldest sample leaving the circular buffer, compared to the running
* average.
*
* Designed to operate with a DC bias inherent in the sample stream.
*
* Tracks a long average and a short average to compensate for drifting
* transmitters.
*
* The long buffer compensates for long strings of ones or
* zeros. Set the long buffer size to a few baud periods longer than the
* longest anticipated string of like symbols.
*
*/
public FloatToBinarySlicer( int bufferSize, int sampleCountForAveraging )
{
mAveragingBuffer = new FloatCircularAveragingBuffer( bufferSize,
sampleCountForAveraging );
}
@Override
public void receive( Float newestSample )
{
float oldestSample = mAveragingBuffer.get( newestSample );
//Dispatch the decision
mListener.receive( oldestSample > mAveragingBuffer.getAverage() );
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void clearListener()
{
mListener = null;
}
}

View File

@@ -0,0 +1,137 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
import sample.Listener;
import buffer.FloatCircularAveragingBuffer2;
/**
* Binary slicer tuned to operate on LTR type short burst signals in a heavy
* DC bias environment coupled with transmitter automatic frequency correction
* being applied as the transmission progresses.
*
* Uses a pair of float buffers that calculate the average of the samples
* contained in the buffer.
*
* The leading and trailing buffers are the same length as the sync pattern.
*
* The slicer uses the Leading buffer's average as the slicing reference for
* normal operation.
*
* Once a sync pattern is detected, the trailing buffer should have the sync
* pattern bauds and the leading buffer should have the first of the message
* bits. At this point, ie sync detect, we switch to using the average of both
* buffers as the slicing reference point.
*
* Once the reference counter counts down to indicate that there are sync length
* bauds remaining of the message, and those remaining bauds should be in the
* leading buffer, we swich to using the trailing buffer as the slicing
* reference point for the remainder of the message.
*
* As the reference counter decrements to 0, meaning the message should be
* clear of both buffers, we switch back to normal operation of using the
* leading buffer as the slicing reference.
*
* An external sync detector indicates a sync pattern match by using the
* syncDetected() method, which immediately switches the slicing average
* reference to both, and sets the counter to the expected message length.
*/
public class FloatToBinarySlicerLTR implements Listener<Float>,
SyncDetectListener
{
private enum SlicingReference
{
LEADING, TRAILING, BOTH;
}
private SlicingReference mSlicingReference = SlicingReference.LEADING;
private FloatCircularAveragingBuffer2 mLeadingBuffer;
private FloatCircularAveragingBuffer2 mTrailingBuffer;
private int mReferenceCounter = 0;
private int mSyncLength;
private int mMessageLength;
private Listener<Boolean> mListener;
public FloatToBinarySlicerLTR( int syncLength,
int messageLength,
int samplesPerBaud )
{
mSyncLength = syncLength * samplesPerBaud;
mMessageLength = messageLength * samplesPerBaud;
mLeadingBuffer = new FloatCircularAveragingBuffer2( mSyncLength );
mTrailingBuffer = new FloatCircularAveragingBuffer2( mSyncLength );
}
@Override
public void receive( Float newestSample )
{
float middleSample = mLeadingBuffer.get( newestSample );
mTrailingBuffer.get( middleSample );
//Dispatch the decision
if( mListener != null )
{
switch( mSlicingReference )
{
case LEADING:
mListener.receive( middleSample > mLeadingBuffer.average() );
break;
case BOTH:
mListener.receive( middleSample >
( ( mLeadingBuffer.average() +
mTrailingBuffer.average() ) / 2.0f ) );
break;
case TRAILING:
mListener.receive( middleSample > mTrailingBuffer.average() );
break;
}
mReferenceCounter--;
if( mReferenceCounter <= mSyncLength && mReferenceCounter > 0 )
{
mSlicingReference = SlicingReference.TRAILING;
}
else if( mReferenceCounter <= 0 )
{
mSlicingReference = SlicingReference.LEADING;
mReferenceCounter = 0;
}
}
}
@Override
public void syncDetected()
{
mSlicingReference = SlicingReference.BOTH;
mReferenceCounter = mMessageLength;
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void clearListener()
{
mListener = null;
}
}

77
src/dsp/NRZDecoder.java Normal file
View File

@@ -0,0 +1,77 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
import sample.Listener;
public class NRZDecoder implements Listener<Boolean>
{
public final static boolean MODE_NORMAL = true;
public final static boolean MODE_INVERTED = false;
private Listener<Boolean> mListener;
private boolean mMode = MODE_NORMAL;
private boolean mPrevious;
/**
* Non Return To Zero (NRZ) and NRZ Inverted (NRZ-I) Decoder.
*
* Use default constructor for NRZ output. Use the alternate constructor
* with the MODE_INVERTED argument for NRZ-I output.
*
* Performs XOR of incoming bit with previous output and sends the result
* to the registered listener.
*/
public NRZDecoder( boolean mode )
{
mMode = mode;
}
public NRZDecoder()
{
this( MODE_NORMAL );
}
public void dispose()
{
mListener = null;
}
@Override
public void receive( Boolean bit )
{
boolean result = mPrevious ^ bit;
if( mListener != null )
{
mListener.receive( mMode ? result : !result );
}
mPrevious = result;
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void removeListener( Listener<Boolean> listener )
{
mListener = null;
}
}

View File

@@ -0,0 +1,27 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
/**
* Interface to allow an external process indicate that it has detected a
* sync pattern in the bit stream
*/
public interface SyncDetectListener
{
public void syncDetected();
}

View File

@@ -0,0 +1,23 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp;
public interface SyncDetectProvider
{
public void setSyncDetectListener( SyncDetectListener listener );
}

View File

@@ -0,0 +1,110 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.Arrays;
import sample.Listener;
/**
* Averages (low pass filter) boolean values over the defined length
*/
public class BooleanAveragingFilter implements Listener<Boolean>
{
private boolean[] mBuffer;
private int mBufferPointer;
private int mThreshold;
private Listener<Boolean> mListener;
public BooleanAveragingFilter( int length )
{
mBuffer = new boolean[ length ];
if( length % 2 == 0 )
{
mThreshold = (int)( length /2 );
}
else
{
mThreshold = (int)( ( length + 1 ) / 2 );
}
//Preload the array with false values
Arrays.fill( mBuffer, false );
}
/**
* Loads the newValue into this buffer and adjusts the buffer pointer
* to prepare for the next get/put cycle
*/
public void put( boolean newValue )
{
//Store the new value to the buffer
mBuffer[ mBufferPointer ] = newValue;
//Increment the buffer pointer
mBufferPointer++;
//Wrap the buffer pointer around to 0 when necessary
if( mBufferPointer >= mBuffer.length )
{
mBufferPointer = 0;
}
}
/**
* Loads the newValue into the buffer, calculates the average
* and returns that average from this method
*
* This effectively performs low-pass filtering
*/
public void receive( Boolean newValue )
{
//Load the new value into the buffer
put( newValue );
int trueCount = 0;
for( int x = 0; x < mBuffer.length; x++ )
{
if( mBuffer[ x ] )
{
trueCount++;
}
}
/**
* If the number of true values in the buffer is
* more than half, send true, otherwise, false
*/
if( mListener != null )
{
mListener.receive( trueCount >= mThreshold );
}
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void clearListener()
{
mListener = null;
}
}

View File

@@ -0,0 +1,341 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
import sample.Provider;
import sample.complex.ComplexSample;
import dsp.filter.Window.WindowType;
/**
* Cascaded Integrator Comb (CIC) filter - anti-aliasing decimating filter for
* reducing a high rate complex sample stream to a lower rate by an integral
* decimation (ie reduction) value.
*
* Note: comb delay is fixed at 1 for this filter.
*/
public class ComplexCICDecimate extends ComplexFilter
implements Provider<ComplexSample>
{
private Integrator[] mIntegrator;
private Comb[] mComb;
private Decimator mDecimator;
private Output mOutput;
/**
* Constructs the complex sample CIC filter
*
* @param decimation - decimation factor
*
* @param stageCount - number of integrator/comb sections to use. The
* higher the number of stages, the sharper dropoff of the filter response
*
* @param outputSampleRate - output sample rate is used to calculate the
* coefficients for the final cleanup filter to adjust for frequency
* response droop induced by this CIC filter
*/
public ComplexCICDecimate( int decimation,
int stageCount,
int outputSampleRate )
{
/* Construct integrator(s) and comb(s) */
mIntegrator = new Integrator[ stageCount ];
mComb = new Comb[ stageCount ];
for( int x = 0; x < stageCount; x++ )
{
mIntegrator[ x ] = new Integrator();
mComb[ x ] = new Comb();
}
/* Connect each stage */
int y = 1;
while( y < stageCount )
{
mIntegrator[ y - 1 ].setListener( mIntegrator[ y ] );
mComb[ y - 1 ].setListener( mComb[ y ] );
y++;
}
/* Connect decimator to integrator output and comb input */
mDecimator = new Decimator( decimation );
mIntegrator[ stageCount - 1 ].setListener( mDecimator );
mDecimator.setListener( mComb[ 0 ] );
/* Connect output adapter to the comb output */
mOutput = new Output( outputSampleRate, stageCount, decimation );
mComb[ stageCount - 1 ].setListener( mOutput );
}
/**
* Sets a new decimation rate for this filter. This method can be invoked
* after filter construction to change the decimation rate dynamically.
*/
public void setDecimationRate( int rate )
{
mDecimator.setDecimationRate( rate );
mOutput.setDecimationRate( rate );
}
/**
* Input method for this filter to receive complex samples.
*/
@Override
public void receive( ComplexSample sample )
{
mIntegrator[ 0 ].receive( (double)sample.left(), (double)sample.right() );
}
/**
* Sets the listener to receive the output of this CIC filter
*/
@Override
public void setListener( Listener<ComplexSample> listener )
{
mOutput.setListener( listener );
}
/**
* Removes the listener from receiving the output of this CIC filter
*/
@Override
public void removeListener( Listener<ComplexSample> listener )
{
mOutput.removeListener( listener );
}
@Override
public Listener<ComplexSample> getListener()
{
return mOutput.getListener();
}
/**
* CIC Comb filter section
*/
public class Comb implements InternalListener
{
private double mPreviousLeft;
private double mPreviousRight;
private InternalListener mCombListener;
/**
* Performs comb filtering on the received values, subtracting the
* previous received samples from these samples.
*/
@Override
public void receive( double left, double right )
{
mCombListener.receive( ( left - mPreviousLeft ),
( right - mPreviousRight ) );
mPreviousLeft = left;
mPreviousRight = right;
}
/**
* Sets the comb output listener.
*/
public void setListener( InternalListener listener )
{
mCombListener = listener;
}
}
/**
* CIC Integrator section
*/
public class Integrator implements InternalListener
{
private double mRunningSumLeft = 0.0d;
private double mRunningSumRight = 0.0d;
private InternalListener mIntegratorListener;
/**
* Integrates the received values with the previous running sum and
* sends the new running sum to the registered listener
*/
@Override
public void receive( double left, double right )
{
mRunningSumLeft += left;
mRunningSumRight += right;
mIntegratorListener.receive( mRunningSumLeft, mRunningSumRight );
}
/**
* Sets the integrator output listener
*/
public void setListener( InternalListener listener )
{
mIntegratorListener = listener;
}
}
/**
* Decimator - passes only 1 of xx samples to a registered listener and
* discards all other samples.
*/
public class Decimator implements InternalListener
{
private InternalListener mDecimatorListener;
private int mPointer = 0;
private int mDecimationRate;
/**
* Constructs a new Decimator object with the specified decimation rate.
*/
public Decimator( int rate )
{
mDecimationRate = rate;
}
/**
* Sets a new decimation rate. Rate can be changed after this object
* is constructed.
*/
public synchronized void setDecimationRate( int rate )
{
mDecimationRate = rate;
}
/**
* Receives samples allowing only 1 of every (rate) sample to go on
* to the registered listener
*/
@Override
public void receive( double left, double right )
{
mPointer++;
if( mPointer >= mDecimationRate )
{
mDecimatorListener.receive( left, right );
mPointer = 0;
}
}
/**
* Sets the decimated output listener
*/
public void setListener( InternalListener listener )
{
mDecimatorListener = listener;
}
}
/**
* Output adapter - applies gain correction, applies cleanup filter, casts
* double-valued samples into float-valued complex sample and sends the
* result to the registered listener.
*/
public class Output implements InternalListener
{
private int mStageCount;
private int mDecimation;
private double mGain = 1.0d;
private ComplexFilter mCleanupFilter;
public Output( int outputSampleRate, int stageCount, int decimation )
{
mStageCount = stageCount;
mDecimation = decimation;
mCleanupFilter = new ComplexFIRFilter( FilterFactory
.getCICCleanupFilter( outputSampleRate,
stageCount,
WindowType.BLACKMAN ), 1.0d );
setGain();
}
/**
* Receiver method for the output adapter to receive a filtered,
* decimated sample, apply gain correction, apply cleanup filtering
* and output the sample values as a complex sample.
*/
@Override
public void receive( double left, double right )
{
mCleanupFilter.receive(
new ComplexSample( (float)( left * mGain ),
(float)( right * mGain ) ) );
}
/**
* Sets a new decimation rate
*/
public void setDecimationRate( int rate )
{
mDecimation = rate;
setGain();
}
/**
* Adds a listener to receive output samples
*/
public void setListener( Listener<ComplexSample> listener )
{
mCleanupFilter.setListener( listener );
}
/**
* Removes the listener from receiving output samples
*/
public void removeListener( Listener<ComplexSample> listener )
{
mCleanupFilter.setListener( null );
}
public Listener<ComplexSample> getListener()
{
return mCleanupFilter.getListener();
}
/**
* Adjusts the gain correction factor for output samples.
*
* The CIC decimation filter induces significant gain as a factor of
* delay, decimation rate and stage count. So, we apply a gain reduction
* to offset the added gain estimate:
*
* gain = ( decimation * delay ) to power of stageCount
*
* Delay is ONE for the Comb used in this filter.
*/
private void setGain()
{
/* decimation is NOT multiplied by delay since delay is 1 */
double gain = Math.pow( (double)mDecimation, (double)mStageCount );
mGain = 10.0d / gain;
}
}
/**
* Internal interface for each of the CIC components to pass samples
*/
private interface InternalListener
{
public void receive( double left, double right );
}
}

View File

@@ -0,0 +1,160 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.complex.ComplexSample;
public class ComplexFIRFilter extends ComplexFilter
{
private ArrayList<ComplexSample> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
private double mGain;
public ComplexFIRFilter( double[] coefficients, double gain )
{
mCoefficients = coefficients;
mBuffer = new ArrayList<ComplexSample>();
mBufferSize = mCoefficients.length;
mGain = gain;
//Fill the buffer with zero valued samples, so we don't have to check for null
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( new ComplexSample( (float)0.0, (float)0.0 ) );
}
generateIndexMap( mCoefficients.length );
}
public void dispose()
{
mBuffer.clear();
mBuffer = null;
mCoefficients = null;
mIndexMap = null;
}
public void receive( ComplexSample newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double leftAccumulator = 0;
double rightAccumulator = 0;
//Start with the center tap value
leftAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).left();
rightAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).right();
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficient; x++ )
{
leftAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).left() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + mCenterCoefficient ] ).left() );
rightAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).right() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + mCenterCoefficient ] ).right() );
}
//We're almost finished ... apply gain, cast the doubles to floats and
//send it on it's merry way
send( new ComplexSample( (float)( leftAccumulator * mGain ),
(float)( rightAccumulator * mGain ) ) );
}
/**
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
mIndexMap = new int[ size ][ size ];
//Last column will be the center coefficient index value for each row
mCenterCoefficientMapIndex = size - 1;
mCenterCoefficient = (int)( size / 2 );
//Setup the first row. Offset is the first row's center coefficient value
//and will become the index value we place in the last column
mIndexMap[ 0 ][ mCenterCoefficientMapIndex ] = mCenterCoefficient;
//Indexes 0 to 1/2 are their index value, and indexes 1/2 to end are
//an offset of the first half of the values
for( int x = 0; x < mCenterCoefficient; x++ )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + mCenterCoefficient ] = size - 1 - x;
}
//For each subsequent map row, increment the value of the preceding row
//same column by 1, wrapping to zero when we exceed the size value
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < size; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] = 0;
}
}
}
}
public String printIndexMap()
{
StringBuilder sb = new StringBuilder();
for( int x = 0; x < mCoefficients.length; x++ )
{
for( int y = 0; y < mCoefficients.length; y++ )
{
sb.append( mIndexMap[ x ][ y ] + " " );
}
sb.append( "\n" );
}
return sb.toString();
}
}

View File

@@ -0,0 +1,54 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
import sample.complex.ComplexSample;
public abstract class ComplexFilter implements Listener<ComplexSample>
{
private Listener<ComplexSample> mListener;
public void setListener( Listener<ComplexSample> listener )
{
mListener = listener;
}
public void dispose()
{
mListener = null;
}
public Listener<ComplexSample> getListener()
{
return mListener;
}
public boolean hasListener()
{
return mListener != null;
}
protected void send( ComplexSample sample )
{
if( mListener != null )
{
mListener.receive( sample );
}
}
}

View File

@@ -0,0 +1,186 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.complex.ComplexSample;
public class ComplexHalfBandFilter extends ComplexFilter
{
private ArrayList<ComplexSample> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double mGain;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
private boolean mDispatchFlag = false;
/**
* Half-Band filter with decimation by 2 against complex samples composed of
* real valued floats.
*
* Takes advantage of the symmetrical nature of FIR filter coefficients by
* adding oldest and newest sample first, then multiplying once by the
* corresponding coefficient
*
* Also, takes advantage of the 0-valued FIR half-band coefficents inherent
* in the half-band filter, and does not calculate those coefficients.
*
* This reduces the workload to (tap-size - 1) / 4 + 1 calculations per sample.
*
* @param filter - filter coefficients
* @param gain - gain multiplier. Use 1.0 for unity/no gain
*/
public ComplexHalfBandFilter( Filters filter, double gain )
{
mCoefficients = filter.getCoefficients();
mBuffer = new ArrayList<ComplexSample>();
mBufferSize = mCoefficients.length;
//Fill the buffer with zero valued samples
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( new ComplexSample( (float)0.0, (float)0.0 ) );
}
generateIndexMap( mCoefficients.length );
mGain = gain;
}
/**
* Calculate the filtered value by applying the coefficients against
* the complex samples in mBuffer
*/
public void receive( ComplexSample newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Toggle the flag every time, so that we can calculate and dispatch a
//sample every other time, when the flag is true (ie decimate by 2)
mDispatchFlag = !mDispatchFlag;
//Calculate a filtered sample when the flag is true
if( mDispatchFlag )
{
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double leftAccumulator = 0;
double rightAccumulator = 0;
//Start with the center tap value
leftAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).left();
rightAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).right();
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficientMapIndex; x += 2 )
{
leftAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).left() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ).left() );
rightAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).right() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ).right() );
}
//We're almost finished ... apply gain, cast the doubles to floats and
//send it on it's merry way
send( new ComplexSample( (float)( leftAccumulator * mGain ),
(float)( rightAccumulator * mGain ) ) );
}
}
/**
* Creates an n X (n + 1 / 2) index map enabling quick access to the
* circular buffer samples.
*
* As the buffer shifts right with each subsequent sample, we have to move
* the index pointers with it, for efficient access of the samples.
*
* The first array index value in the index map corresponds to the current
* buffer pointer location.
*
* The second array index value points to the samples that should be
* multiplied by the coefficients as follows:
*
* 0 = center tap sample, to be multiplied by center coefficient
*
* 0 = sample( 1 )
* 1 = sample( size - 1 )
*
* Indexes 0 and 1 will be multiplied by coefficient( 0 ).
*
* Subsequent indexes 3, 4, etc, point to the oldest and newest samples that
* correspond to the matching ( 3 ) coefficient index.
*
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
int mapWidth = ( ( size + 1 ) / 2 ) + 1;
//Set center tap index for coefficients array
mCenterCoefficient = ( size - 1 ) / 2;
mCenterCoefficientMapIndex = mCenterCoefficient + 1;
mIndexMap = new int[ size ][ mapWidth ];
//Setup the first row, buffer pointer index 0, as a starting point
for( int x = 0; x < mapWidth - 2; x += 2 )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + 1 ] = size - 1 - x;
}
//Place center tap index in last element
mIndexMap[ 0 ][ mCenterCoefficientMapIndex ] = mCenterCoefficient;
//For each subsequent row, increment the previous row's value by 1,
//subtracting size as needed, to keep the values between 0 and size - 1
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < mapWidth; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] -= size;
}
}
}
}
}

View File

@@ -0,0 +1,180 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.Listener;
import sample.complex.ComplexSample;
public class ComplexHalfBandNoDecimateFilter extends ComplexFilter
{
private Listener<ComplexSample> mListener;
private ArrayList<ComplexSample> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double mGain;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
/**
* Half-Band filter with no decimation against complex samples composed of
* real valued floats.
*
* Takes advantage of the symmetrical nature of FIR filter coefficients by
* adding oldest and newest sample first, then multiplying once by the
* corresponding coefficient
*
* Also, takes advantage of the 0-valued FIR half-band coefficents inherent
* in the half-band filter, and does not calculate those coefficients.
*
* This reduces the workload to (tap-size - 1) / 4 + 1 calculations per sample.
*
* @param filter - filter coefficients
* @param gain - gain multiplier. Use 1.0 for unity/no gain
*/
public ComplexHalfBandNoDecimateFilter( Filters filter, double gain )
{
mCoefficients = filter.getCoefficients();
mBuffer = new ArrayList<ComplexSample>();
mBufferSize = mCoefficients.length;
//Fill the buffer with zero valued samples
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( new ComplexSample( (float)0.0, (float)0.0 ) );
}
generateIndexMap( mCoefficients.length );
mGain = gain;
}
/**
* Calculate the filtered value by applying the coefficients against
* the complex samples in mBuffer
*/
public void receive( ComplexSample newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double leftAccumulator = 0;
double rightAccumulator = 0;
//Start with the center tap value
leftAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).left();
rightAccumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] ).right();
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficientMapIndex; x += 2 )
{
leftAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).left() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ).left() );
rightAccumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ).right() +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ).right() );
}
//We're almost finished ... apply gain, cast the doubles to floats and
//send it on it's merry way
send( new ComplexSample( (float)( leftAccumulator * mGain ),
(float)( rightAccumulator * mGain ) ) );
}
/**
* Creates an n X (n + 1 / 2) index map enabling quick access to the
* circular buffer samples.
*
* As the buffer shifts right with each subsequent sample, we have to move
* the index pointers with it, for efficient access of the samples.
*
* The first array index value in the index map corresponds to the current
* buffer pointer location.
*
* The second array index value points to the samples that should be
* multiplied by the coefficients as follows:
*
* 0 = center tap sample, to be multiplied by center coefficient
*
* 0 = sample( 1 )
* 1 = sample( size - 1 )
*
* Indexes 0 and 1 will be multiplied by coefficient( 0 ).
*
* Subsequent indexes 3, 4, etc, point to the oldest and newest samples that
* correspond to the matching ( 3 ) coefficient index.
*
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
int mapWidth = ( ( size + 1 ) / 2 ) + 1;
//Set center tap index for coefficients array
mCenterCoefficient = ( size - 1 ) / 2;
mCenterCoefficientMapIndex = mCenterCoefficient + 1;
mIndexMap = new int[ size ][ mapWidth ];
//Setup the first row, buffer pointer index 0, as a starting point
for( int x = 0; x < mapWidth - 2; x += 2 )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + 1 ] = size - 1 - x;
}
//Place center tap index in last element
mIndexMap[ 0 ][ mCenterCoefficientMapIndex ] = mCenterCoefficient;
//For each subsequent row, increment the previous row's value by 1,
//subtracting size as needed, to keep the values between 0 and size - 1
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < mapWidth; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] -= size;
}
}
}
}
}

View File

@@ -0,0 +1,177 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
import sample.complex.ComplexSample;
import buffer.ComplexTappedCircularBuffer;
public class ComplexPolyphaseDecimationFilter extends ComplexFilter
{
private FilterPhase[] mFilterPhases;
private int mPhasePointer = 0;
private int mPhaseCount;
/**
* Implements a decimating poly-phase filter bank for complex samples, as
* described in Understanding Digital Signal Processing, 3e, Lyons, in
* section 10.7, figure 10-14.
*
* Constraints:
* - Low-pass filter coefficient length must be an integer multiple
* of the decimation rate:
*
* Consideration:
* - Low pass filter cutoff frequency must be less than:
* ( input sample rate / 2 / decimation )
*
* @param coefficients - low pass filter coefficients.
* @param decimationRate - integer decimation rate
*/
public ComplexPolyphaseDecimationFilter( double[] coefficients,
int decimationRate )
{
/* Verify that the filter length is an integer multiple of the
* decimation rate */
assert( coefficients.length % decimationRate == 0 );
mPhaseCount = decimationRate;
/* Construct a number of filter phases equal to the decimation rate */
mFilterPhases = new FilterPhase[ mPhaseCount ];
for( int phaseNumber = 0; phaseNumber < mPhaseCount; phaseNumber++ )
{
double[] phaseCoefficients =
getPhaseCoefficients( coefficients, decimationRate, phaseNumber );
mFilterPhases[ phaseNumber ] = new FilterPhase( phaseCoefficients );
}
}
/**
* Inputs a new sample into the filter at the non-decimated rate
*/
@Override
public void receive( ComplexSample sample )
{
mFilterPhases[ mPhasePointer++ ].receive( sample );
if( mPhasePointer >= mPhaseCount )
{
/* Broadcast a new output sample */
calculate();
/* Reset the pointer */
mPhasePointer = 0;
}
}
/**
* Outputs a complex sample at the decimated rate, comprised of the sum of
* the outputs of the individual filter stages, trigger by a decimation rate
* counter on the input side.
*/
private void calculate()
{
if( hasListener() )
{
/* Accumulators */
float left = 0.0f;
float right = 0.0f;
for( int x = 0; x < mPhaseCount; x++ )
{
ComplexSample phaseSample = mFilterPhases[ x ].calculate();
left += phaseSample.left();
right += phaseSample.right();
}
send( new ComplexSample( left, right ) );
}
}
/**
* Returns the subset of filter coefficients that belong to the filter bank
* phase.
*/
private double[] getPhaseCoefficients( double[] coefficients,
int phaseCount,
int phaseNumber )
{
int phaseFilterLength = (int)( coefficients.length / phaseCount );
double[] phaseCoefficients = new double[ phaseFilterLength ];
for( int x = 0; x < phaseFilterLength; x++ )
{
phaseCoefficients[ x ] = coefficients[ x * phaseCount + phaseNumber ];
}
return phaseCoefficients;
}
/**
* Implements a single phase of a poly-phase filter bank
*/
public class FilterPhase implements Listener<ComplexSample>
{
private int mLength;
private double[] mCoefficients;
private ComplexTappedCircularBuffer mBuffer;
/**
* @param coefficients - subset of the filter coefficients for this
* phase of the poly-phase filter bank
*/
public FilterPhase( double[] coefficients )
{
mCoefficients = coefficients;
mLength = coefficients.length;
/* The constructor preloads the buffer with 0-valued samples */
mBuffer = new ComplexTappedCircularBuffer( coefficients.length );
}
@Override
public void receive( ComplexSample sample )
{
mBuffer.add( sample );
}
/**
* Performs convolution by multiplying the buffer samples by their
* corresponding filter coefficients and returns the sum
*/
public ComplexSample calculate()
{
/* Accumulators */
double left = 0.0d;
double right = 0.0d;
for( int x = 0; x < mLength; x++ )
{
left += ( mBuffer.get( x ).left() * mCoefficients[ x ] );
right += ( mBuffer.get( x ).right() * mCoefficients[ x ] );
}
return new ComplexSample( (float)left, (float)right );
}
}
}

View File

@@ -0,0 +1,65 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
/**
* Implements the real-time general DC removal algorithm described on page
* 762, figure 13-62d, in Understanding Digital Signal Processing 3e, Lyons.
*/
public class DCRemovalFilter implements Listener<Float>
{
private double mAlpha;
private float mPrevious;
private Listener<Float> mListener;
public DCRemovalFilter( double alpha )
{
mAlpha = alpha;
}
public void dispose()
{
mListener = null;
}
@Override
public void receive( Float sample )
{
send( (float)( sample + ( mPrevious * mAlpha ) - mPrevious ) );
mPrevious = sample;
}
/**
* Sends the filtered sample to the listener
*/
private void send( float sample )
{
if( mListener != null )
{
mListener.receive( sample );
}
}
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
}

View File

@@ -0,0 +1,56 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
public class DCRemovalFilter2 implements Listener<Float>
{
private double mAverage;
private double mRatio;
private Listener<Float> mListener;
public DCRemovalFilter2( double ratio )
{
mRatio = ratio;
}
@Override
public void receive( Float sample )
{
mAverage += mRatio * ( (double)sample - mAverage );
send( (float)( sample - mAverage ) );
}
/**
* Sends the filtered sample to all registered listeners
*/
private void send( float sample )
{
if( mListener != null )
{
mListener.receive( sample );
}
}
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
}

View File

@@ -0,0 +1,729 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.Arrays;
import java.util.Set;
import java.util.TreeSet;
import log.Log;
import dsp.filter.Window.WindowType;
import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
public class FilterFactory
{
public static void main( String[] args )
{
Log.info( "Testing filter factory" );
double[] coefficients = getCICCleanupFilter( 48000, 1, WindowType.BLACKMAN );
for( int x = 0; x < coefficients.length; x++ )
{
Log.info( x + ": " + coefficients[ x ] );
}
Log.info( "Finished" );
}
/**
* Generates coefficients for a unity-gain, windowed low-pass filter
* @param sampleRate - hertz
* @param frequency - cutoff in hertz
* @param length - filter length
* @param window - to apply against the coefficients
* @return
*/
public static double[] getSinc( int sampleRate,
int frequency,
int length,
WindowType window )
{
//Ensure we have an odd length
assert( length % 2 == 0 );
//Get unity response array (one element longer to align with IDFT size)
double[] frequencyResponse = getUnityResponseArray( sampleRate, frequency, length + 1 );
//Apply Inverse DFT against frequency response unity values, leaving the
//IDFT bin results in the frequency response array
DoubleFFT_1D idft = new DoubleFFT_1D( length + 1 );
idft.realInverseFull( frequencyResponse, false );
//Transfer the IDFT results to the return array
double[] coefficients = new double[ length ];
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++ )
{
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;
}
/**
* Normalizes all filter coefficients to achieve unity (1) gain, by ensuring
* that the sum of the absolute value of all coefficients adds up to 1.
* @param coefficients
* @return
*/
private static double[] normalize( double[] coefficients )
{
double accumulator = 0;
for( int x = 0; x < coefficients.length; x++ )
{
accumulator += Math.abs( coefficients[ x ] );
}
for( int x = 0; x < coefficients.length; x++ )
{
coefficients[ x ] = coefficients[ x ] / accumulator;
}
return coefficients;
}
/**
* Constructs an array of unity (1) response values representing the
* desired (pre-windowing) frequency response, used by the sync function
*
* Returns an array twice the length, filled with unity (1) response values
* in the desired pass-band with the positive frequency response starting at
* the lower end of the array, and the negative frequency response at the
* higher end of the first half of the array. The remaining zero-valued
* indexes in the second half will store the results of the JTransforms
* inverse DFT operation
*
* @param sampleRate
* @param frequency
* @param length
* @return
*/
public static double[] getUnityResponseArray( int sampleRate,
int frequency,
int length )
{
double[] unityArray = new double[ length * 2 ];
int binCount = (int)( ( Math.round(
(double)frequency / (double)sampleRate * (double)length ) ) );
// if( length % 2 == 0 ) //even length
// {
// for( int x = 0; x < binCount; x++ )
// {
// unityArray[ x ] = 1.0d;
// unityArray[ length - 1 - x ] = 1.0d;
// }
// }
// else //odd length
// {
unityArray[ 0 ] = 1.0d;
for( int x = 1; x <= binCount; x++ )
{
unityArray[ x ] = 1.0d;
unityArray[ length - x ] = 1.0d;
}
// }
return unityArray;
}
public static double[] getSine( double sampleRate, double frequency, int length )
{
//Ensure we have an odd length
assert( length % 2 == 0 );
double[] retVal = new double[ length ];
double radianFrequency = 2.0D * Math.PI * ( frequency / sampleRate );
int middle = (int)( length / 2 );
for( int x = 0; x < middle; x++ )
{
double val = 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
* filter to a high-pass filter.
*
* Inverts the sign of all odd index coefficients ( 1, 3, 5, etc.)
* returning:
* Index 0: same
* Index 1: inverted
* Index 2: same
* ...
* Index length - 1: same
*/
public static double[] invert( double[] coefficients )
{
for( int x = 1; x < coefficients.length; x += 2 )
{
coefficients[ x ] = -coefficients[ x ];
}
return coefficients;
}
/**
* Generates filter coefficients for a unity-gain, odd-length, windowed,
* low pass filter with passband from 0-hertz to the cutoff frequency.
*
* @param sampleRate - hertz
* @param cutoff - frequency in hertz
* @param filterLength - odd filter length
* @param windowType - window to apply against the generated coefficients
* @return
*/
public static double[] getLowPass( int sampleRate,
int cutoff,
int filterLength,
WindowType windowType )
{
if( filterLength % 2 == 0 ) //even length
{
double[] 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
{
double[] values = getSinc( sampleRate, cutoff, filterLength + 1, windowType );
//throw away the 0 index
return Arrays.copyOfRange( values, 1, values.length );
}
}
/**
* Creates a low-pass filter with ~ 0.1 db ripple in the pass band
*
* Note: stop frequency - pass frequency defines the transition band.
*/
/**
* Creates a low-pass filter with ~0.1 db ripple in the pass band. The
* transition region is defined by the stop frequency minus the pass
* frequency.
*
* Requires:
* - passFrequency < stopFrequency
* - stopFrequency <= sampleRate/2
*/
public static double[] getLowPass( int sampleRate,
int passFrequency,
int stopFrequency,
int attenuation,
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 + "]" );
}
int tapCount = getTapCount( sampleRate,
passFrequency,
stopFrequency,
attenuation );
if( forceOddLength )
{
if( tapCount % 2 == 0 )
{
tapCount--;
}
}
return getLowPass( sampleRate, passFrequency, tapCount, windowType );
}
/**
* Generates filter coefficients for a unity-gain, odd-length, windowed,
* high pass filter with passband from cutoff frequency to half the sample
* rate.
*
* @param sampleRate - hertz
* @param cutoff - frequency in hertz
* @param filterLength - odd filter length
* @param windowType - window to apply against the generated coefficients
* @return
*/
public static double[] getHighPass( int sampleRate,
int cutoff,
int filterLength,
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
int convertedCutoff = sampleRate / 2 - cutoff;
return invert( getSinc( sampleRate,
convertedCutoff,
filterLength,
windowType ) );
}
public static double[] getHighPass( int sampleRate,
int stopFrequency,
int passFrequency,
int attenuation,
WindowType windowType,
boolean forceOddLength )
{
/* reverse the stop and pass frequency to get the low pass variant */
int tapCount = getTapCount( sampleRate,
stopFrequency,
passFrequency,
attenuation );
if( forceOddLength )
{
if( tapCount % 2 == 0 )
{
tapCount--;
}
}
return invert( getLowPass( sampleRate, stopFrequency, tapCount, windowType ) );
}
/**
* Utility to log the arrays of doubles with line breaks
*/
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 ] );
if( breaks )
{
sb.append( "\n" );
}
else
{
if( x % 8 == 7 )
{
sb.append( "\n" );
}
else
{
sb.append( "\t" );
}
}
}
return sb.toString();
}
/**
* Determines the number of fir filter taps required to produce the
* specified frequency response with passband ripple near .1dB.
*
* Implements the algorithm from Understanding Digital Signal Processing, 3e
* , Lyons, section 5.10.5.
*
* @param sampleRate in hertz
* @param pass pass frequency in hertz
* @param stop stop frequency in hertz
* @param attenuation in dB
* @return
*/
public static int getTapCount( int sampleRate,
int pass,
int stop,
int attenuation )
{
double frequency = ( (double)stop - (double)pass ) / (double)sampleRate;
return (int)( Math.round( (double)attenuation / ( 22.0d * frequency ) ) );
}
public static ComplexFilter[] getDecimationFilters( int sampleRate,
int decimatedRate, int passFrequency, int attenuation,
WindowType windowType )
{
ComplexFilter[] filters;
if( sampleRate == 96000 && decimatedRate == 48000 )
{
filters = new ComplexFilter[ 1 ];
filters[ 0 ] = new ComplexHalfBandFilter(
Filters.FIR_HALF_BAND_31T_ONE_EIGHTH_FCO, 1.0002 );
}
else if( sampleRate == 192000 && decimatedRate == 48000 )
{
filters = new ComplexFilter[ 2 ];
filters[ 0 ] = new ComplexHalfBandFilter(
Filters.FIR_HALF_BAND_31T_ONE_EIGHTH_FCO, 1.0002 );
filters[ 1 ] = new ComplexHalfBandFilter(
Filters.FIR_HALF_BAND_31T_ONE_EIGHTH_FCO, 1.0002 );
}
else
{
filters = new ComplexFilter[ 1 ];
int decimate = (int)( sampleRate / decimatedRate );
filters[ 0 ] = new ComplexCICDecimate( decimate, 2, decimatedRate );
// int finalStopFrequency = decimatedRate - passFrequency;
//
// int[] rates = getPolyphaseDecimationRates( sampleRate,
// decimatedRate,
// passFrequency,
// finalStopFrequency );
//
// filters = new ComplexFilter[ rates.length ];
//
// int counter = 0;
// int currentSampleRate = sampleRate;
// int currentDecimatedRate = sampleRate;
//
// while( counter < rates.length )
// {
// currentDecimatedRate /= rates[ counter ];
//
// int stopFrequency = currentDecimatedRate- passFrequency;
//
// int proposedFilterLength = getTapCount( currentSampleRate,
// passFrequency,
// stopFrequency,
// attenuation );
//
// /* check filter length is multiple of decimation rate */
// int multiplier = (int)( proposedFilterLength / rates[ counter ] );
//
// /* increase filter length to the next multiple, if needed */
// if( multiplier * rates[ counter ] < proposedFilterLength )
// {
// multiplier++;
// }
//
// double[] filter = getLowPass( currentSampleRate,
// passFrequency,
// rates[counter] * multiplier,
// windowType );
//
// filters[ counter ] =
// new ComplexPolyphaseDecimationFilter( filter, rates[ counter ] );
//
// Log.info( "FilterFactory - PFD filter [" + ( counter + 1 ) + "/" + rates.length +
// "] overall rate [" + sampleRate +
// "] decimate [" + rates[ counter ] +
// "] filter taps [" + filter.length +
// "] rate [" + currentSampleRate +
// "] pass [" + passFrequency +
// "]" );
//
// currentSampleRate /= rates[ counter ];
// counter++;
// }
}
return filters;
}
/**
* Determines decimation rate(s) for a polyphase decimation filter.
*
* @param sampleRate - starting sample rate
* @param decimatedRate - final (decimated) output rate
* @return - set of integer decimation rates for a single or multi-stage
* polyphase filter decimation chain
* @throws - AssertionException if sample rate is not a multiple of 48 kHz
*/
public static int[] getPolyphaseDecimationRates( int sampleRate,
int decimatedRate,
int passFrequency,
int stopFrequency )
{
int[] rates;
if( sampleRate % decimatedRate != 0 )
{
throw new IllegalArgumentException( "Decimated rate must be an "
+ "integer multiple of sample rate" );
}
int decimation = (int)( sampleRate / decimatedRate );
//Decimation rates below 20 will use single stage polyphase filter
if( decimation < 20 )
{
rates = new int[ 1 ];
rates[ 0 ] = decimation;
return rates;
}
else
{
int optimalStage1 =
getOptimalStageOneRate( sampleRate, decimation, passFrequency, stopFrequency );
Set<Integer> factors = getFactors( decimation );
int stage1 = findClosest( optimalStage1, factors );
Log.info( "Decimation rate [" + decimation +
"] stage1 optimal [" + optimalStage1 +
"] stage1 actual [" + stage1 +
"]");
if( stage1 == decimation || stage1 == 1 )
{
rates = new int[ 1 ];
rates[ 0 ] = decimation;
return rates;
}
else
{
rates = new int[ 2 ];
int stage2 = (int)( decimation / stage1);
if( stage1 > stage2 )
{
rates[ 0 ] = stage1;
rates[ 1 ] = stage2;
}
else
{
rates[ 0 ] = stage2;
rates[ 1 ] = stage1;
}
return rates;
}
}
}
/**
* 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 )
{
int bestFactor = 1;
int bestDelta = desiredFactor;
for( Integer factor: factors )
{
int testDelta = Math.abs( desiredFactor - factor );
if( testDelta < bestDelta )
{
bestDelta = testDelta;
bestFactor = factor;
}
}
return bestFactor;
}
/**
* 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
*/
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 );
if( remainder * x == value )
{
factors.add( x );
}
}
return factors;
}
/**
* 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.
*
* Implements the algorithm described in Lyons, Understanding Digital Signal
* Processing, 3e, section 10.2.1, page 511.
*
* @param sampleRate
* @param decimation
* @param passFrequency frequency of the pass band
* @return optimum integer decimation rate for the first stage decimation
* filter
*/
public static int getOptimalStageOneRate( int sampleRate,
int decimation,
int passFrequency,
int stopFrequency )
{
double ratio = getBandwidthRatio( passFrequency, stopFrequency );
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 ) );
Log.info( "Optimal Stage 1 Decimation - rate [" + sampleRate +
"] pass [" + passFrequency +
"] bw ratio [" + ratio +
"] optimal [" + retVal +
"]" );
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.
*/
private static double getBandwidthRatio( int passFrequency, int stopFrequency )
{
assert( passFrequency < stopFrequency );
return ( (double)( stopFrequency - passFrequency ) /
(double)stopFrequency );
}
/**
* Assumes that the pass band is 1/4 of the output sample rate.
*
* Assumes the stop band is: pass + (pass * .25).
* @param outputSampleRate
* @param stageCount
* @return
*/
public static double[] getCICCleanupFilter( int outputSampleRate,
int stageCount,
WindowType window )
{
int pass = (int)( outputSampleRate / 4 );
int stop = pass + (int)( (double)pass * 0.25d );
int attenuation = 60; //db
int taps = getTapCount( outputSampleRate, pass, stop, attenuation );
/* Make tap count odd */
if( taps% 2 == 0 )
{
taps ++;
}
double[] frequencyResponse =
getCICResponseArray( outputSampleRate, pass, taps, stageCount );
//Apply Inverse DFT against frequency response unity values, leaving the
//IDFT bin results in the frequency response array
DoubleFFT_1D idft = new DoubleFFT_1D( taps );
idft.realInverseFull( frequencyResponse, false );
//Transfer the IDFT results to the odd length return array
double[] coefficients = new double[ taps ];
int middleCoefficient = (int)( taps / 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++ )
{
coefficients[ middleCoefficient + x ] = frequencyResponse[ 2 * x ];
coefficients[ middleCoefficient - x ] = frequencyResponse[ 2 * x ];
}
//Apply the window against the coefficients
coefficients = Window.apply( window, coefficients );
return coefficients;
}
public static double[] getCICResponseArray( int sampleRate,
int frequency,
int length,
int stageCount )
{
double[] cicArray = new double[ length * 2 ];
int binCount = (int)( ( Math.round(
(double)frequency / (double)sampleRate * (double)length ) ) );
cicArray[ 0 ] = 1.0d;
double unityResponse = Math.pow( Math.sin( 1.0d / (double)length ) /
( 1.0d / (double)length ), stageCount );
for( int x = 1; x <= binCount; x++ )
{
/* Calculate unity response plus amplification for droop */
double compensated = 1.0d + ( unityResponse -
Math.pow( ( Math.sin( (double)x / (double)length ) /
( (double)x / (double)length ) ), stageCount ) );
cicArray[ x ] = compensated;
cicArray[ length - x ] = compensated;
}
return cicArray;
}
}

1238
src/dsp/filter/Filters.java Normal file

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,177 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.Listener;
public class FloatFIRFilter implements Listener<Float>
{
private Listener<Float> mListener;
private ArrayList<Float> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
private double mGain;
public FloatFIRFilter( double[] coefficients, double gain )
{
mCoefficients = coefficients;
mBuffer = new ArrayList<Float>();
mBufferSize = mCoefficients.length;
mGain = gain;
//Fill the buffer with zero valued samples, so we don't have to check for null
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( 0.0f );
}
generateIndexMap( mCoefficients.length );
}
public void dispose()
{
mListener = null;
mBuffer.clear();
}
public int getTapCount()
{
return mCoefficients.length;
}
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
public void receive( Float newSample )
{
send( get( newSample ) );
}
public float get( float newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double accumulator = 0.0d;
//Start with the center tap value
accumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] );
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficient; x++ )
{
accumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ) +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + mCenterCoefficient ] ) );
}
//We're almost finished ... apply gain, cast the doubles to shorts and
//return the value
return (float)( accumulator * mGain );
}
/**
* Dispatch the filtered sample to all registered listeners
*/
private void send( float sample )
{
if( mListener != null )
{
mListener.receive( sample );
}
}
public String printIndexMap()
{
StringBuilder sb = new StringBuilder();
for( int x = 0; x < mCoefficients.length; x++ )
{
for( int y = 0; y < mCoefficients.length; y++ )
{
sb.append( mIndexMap[ x ][ y ] + " " );
}
sb.append( "\n" );
}
return sb.toString();
}
/**
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
mIndexMap = new int[ size ][ size ];
//Last column will be the center coefficient index value for each row
mCenterCoefficientMapIndex = size - 1;
mCenterCoefficient = (int)( size / 2 );
//Setup the first row. Offset is the first row's center coefficient value
//and will become the index value we place in the last column
mIndexMap[ 0 ][ mCenterCoefficient ] = mCenterCoefficient;
//Indexes 0 to 1/2 are their index value, and indexes 1/2 to end are
//an offset of the first half of the values
for( int x = 0; x < mCenterCoefficient; x++ )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + mCenterCoefficient ] = size - 1 - x;
}
//For each subsequent map row, increment the value of the preceding row
//same column by 1, wrapping to zero when we exceed the size value
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < size; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] = 0;
}
}
}
}
}

View File

@@ -0,0 +1,203 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.Listener;
public class FloatHalfBandFilter implements Listener<Float>
{
private Listener<Float> mListener;
private ArrayList<Float> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double mGain;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
private boolean mDispatchFlag = false;
/**
* Half-Band filter with decimation by 2 against real valued floats.
*
* Takes advantage of the symmetrical nature of FIR filter coefficients by
* adding oldest and newest sample first, then multiplying once by the
* corresponding coefficient
*
* Also, takes advantage of the 0-valued FIR half-band coefficents inherent
* in the half-band filter, and does not calculate those coefficients.
*
* This reduces the workload to (tap-size - 1) / 4 + 1 calculations per sample.
*
* @param filter - filter coefficients
* @param gain - gain multiplier. Use 1.0 for unity/no gain
*/
public FloatHalfBandFilter( Filters filter, double gain )
{
mCoefficients = filter.getCoefficients();
mBuffer = new ArrayList<Float>();
mBufferSize = mCoefficients.length;
//Fill the buffer with zero valued samples
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( 0.0f );
}
generateIndexMap( mCoefficients.length );
mGain = gain;
}
public void dispose()
{
mListener = null;
}
/**
* Calculate the filtered value by applying the coefficients against
* the complex samples in mBuffer
*/
public void receive( Float newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Toggle the flag every time, so that we can calculate and dispatch a
//sample every other time, when the flag is true (ie decimate by 2)
mDispatchFlag = !mDispatchFlag;
//Calculate a filtered sample when the flag is true
if( mDispatchFlag )
{
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double accumulator = 0;
//Start with the center tap value
accumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] );
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficientMapIndex; x += 2 )
{
accumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ) +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ) );
}
//We're almost finished ... apply gain, cast the doubles to shorts and
//send it on it's merry way
if( mListener != null )
{
mListener.receive( (float)( accumulator * mGain ) );
}
}
}
/**
* Creates an n X (n + 1 / 2) index map enabling quick access to the
* circular buffer samples.
*
* As the buffer shifts right with each subsequent sample, we have to move
* the index pointers with it, for efficient access of the samples.
*
* The first array index value in the index map corresponds to the current
* buffer pointer location.
*
* The second array index value points to the samples that should be
* multiplied by the coefficients as follows:
*
* 0 = center tap sample, to be multiplied by center coefficient
*
* 0 = sample( 1 )
* 1 = sample( size - 1 )
*
* Indexes 0 and 1 will be multiplied by coefficient( 0 ).
*
* Subsequent indexes 3, 4, etc, point to the oldest and newest samples that
* correspond to the matching ( 3 ) coefficient index.
*
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
int mapWidth = ( ( size + 1 ) / 2 ) + 1;
//Set center tap index for coefficients array
mCenterCoefficient = ( size - 1 ) / 2;
mCenterCoefficientMapIndex = mCenterCoefficient + 1;
mIndexMap = new int[ size ][ mapWidth ];
//Setup the first row, buffer pointer index 0, as a starting point
for( int x = 0; x < mapWidth - 2; x += 2 )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + 1 ] = size - 1 - x;
}
//Place center tap index in last element
mIndexMap[ 0 ][ mCenterCoefficientMapIndex ] = mCenterCoefficient;
//For each subsequent row, increment the previous row's value by 1,
//subtracting size as needed, to keep the values between 0 and size - 1
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < mapWidth; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] -= size;
}
}
}
}
/**
* Registers a listener for filtered samples
*/
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
/**
* Removes (if exists) a registered filtered sample listener
*/
public void clearListener()
{
mListener = null;
}
}

View File

@@ -0,0 +1,193 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import java.util.ArrayList;
import sample.Listener;
public class FloatHalfBandNoDecimateFilter implements Listener<Float>
{
private Listener<Float> mListener;
private ArrayList<Float> mBuffer;
private int mBufferSize = 1; //Temporary initial value
private int mBufferPointer = 0;
private double mGain;
private double[] mCoefficients;
private int[][] mIndexMap;
private int mCenterCoefficient;
private int mCenterCoefficientMapIndex;
/**
* Half-Band filter with decimation by a real valued float.
*
* Takes advantage of the symmetrical nature of FIR filter coefficients by
* adding oldest and newest sample first, then multiplying once by the
* corresponding coefficient
*
* Also, takes advantage of the 0-valued FIR half-band coefficents inherent
* in the half-band filter, and does not calculate those coefficients.
*
* This reduces the workload to (tap-size - 1) / 4 + 1 calculations per sample.
*
* @param filter - filter coefficients
* @param gain - gain multiplier. Use 1.0 for unity/no gain
*/
public FloatHalfBandNoDecimateFilter( Filters filter, double gain )
{
mCoefficients = filter.getCoefficients();
mBuffer = new ArrayList<Float>();
mBufferSize = mCoefficients.length;
//Fill the buffer with zero valued samples
for( int x = 0; x < mCoefficients.length; x++ )
{
mBuffer.add( 0.0f );
}
generateIndexMap( mCoefficients.length );
mGain = gain;
}
public void dispose()
{
mListener = null;
}
/**
* Calculate the filtered value by applying the coefficients against
* the complex samples in mBuffer
*/
public void receive( Float newSample )
{
//Add the new sample to the buffer
mBuffer.set( mBufferPointer, newSample );
//Increment & Adjust the buffer pointer for circular wrap around
mBufferPointer++;
if( mBufferPointer >= mBufferSize )
{
mBufferPointer = 0;
}
//Convolution - multiply filter coefficients by the circular buffer
//samples to calculate a new filtered value
double accumulator = 0;
//Start with the center tap value
accumulator += mCoefficients[ mCenterCoefficient ] *
mBuffer.get( mIndexMap[ mBufferPointer ][ mCenterCoefficientMapIndex ] );
//For the remaining coefficients, add the symmetric samples, oldest and newest
//first, then multiply by the single coefficient
for( int x = 0; x < mCenterCoefficientMapIndex; x += 2 )
{
accumulator += mCoefficients[ x ] *
( mBuffer.get( mIndexMap[ mBufferPointer ][ x ] ) +
mBuffer.get( mIndexMap[ mBufferPointer ][ x + 1 ] ) );
}
//We're almost finished ... apply gain, cast the doubles to floats and
//send it on it's merry way
if( mListener != null )
{
mListener.receive( (float)( accumulator * mGain ) );
}
}
/**
* Creates an n X (n + 1 / 2) index map enabling quick access to the
* circular buffer samples.
*
* As the buffer shifts right with each subsequent sample, we have to move
* the index pointers with it, for efficient access of the samples.
*
* The first array index value in the index map corresponds to the current
* buffer pointer location.
*
* The second array index value points to the samples that should be
* multiplied by the coefficients as follows:
*
* 0 = center tap sample, to be multiplied by center coefficient
*
* 0 = sample( 1 )
* 1 = sample( size - 1 )
*
* Indexes 0 and 1 will be multiplied by coefficient( 0 ).
*
* Subsequent indexes 3, 4, etc, point to the oldest and newest samples that
* correspond to the matching ( 3 ) coefficient index.
*
* @param odd-sized number of filter taps (ie coefficients) and buffer
*/
private void generateIndexMap( int size )
{
//Ensure we have an odd size
assert( size % 2 == 1 );
int mapWidth = ( ( size + 1 ) / 2 ) + 1;
//Set center tap index for coefficients array
mCenterCoefficient = ( size - 1 ) / 2;
mCenterCoefficientMapIndex = mCenterCoefficient + 1;
mIndexMap = new int[ size ][ mapWidth ];
//Setup the first row, buffer pointer index 0, as a starting point
for( int x = 0; x < mapWidth - 2; x += 2 )
{
mIndexMap[ 0 ][ x ] = x;
mIndexMap[ 0 ][ x + 1 ] = size - 1 - x;
}
//Place center tap index in last element
mIndexMap[ 0 ][ mCenterCoefficientMapIndex ] = mCenterCoefficient;
//For each subsequent row, increment the previous row's value by 1,
//subtracting size as needed, to keep the values between 0 and size - 1
for( int x = 1; x < size; x++ )
{
for( int y = 0; y < mapWidth; y++ )
{
mIndexMap[ x ][ y ] = mIndexMap[ x - 1 ][ y ] + 1;
if( mIndexMap[ x ][ y ] >= size )
{
mIndexMap[ x ][ y ] -= size;
}
}
}
}
/**
* Registers a listener for filtered samples
*/
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
/**
* Removes (if exists) a registered filtered sample listener
*/
public void clearListener()
{
mListener = null;
}
}

45
src/dsp/filter/Gain.java Normal file
View File

@@ -0,0 +1,45 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
public class Gain implements Listener<Float>
{
private double mGain;
private Listener<Float> mListener;
public Gain( double gain )
{
mGain = gain;
}
@Override
public void receive( Float sample )
{
if( mListener != null )
{
mListener.receive( (float)( sample * mGain ) );
}
}
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
}

View File

@@ -0,0 +1,161 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import dsp.filter.Window.WindowType;
/**
* Implements the Goertzel Filter as described in
* http://en.wikipedia.org/wiki/Goertzel_algorithm
*
* Calculates the magnitude of a specific frequency component within a window of
* samples, using a Goertzel Filter (ie an optimized DFT).
*/
public class GoertzelFilter
{
private int mSampleRate;
private int mTargetFrequency;
private int mBlockSize;
private double[] mWindowCoefficients;
private double mCoefficient;
/**
* Constructs the Goertzel Filter to use for block processing of
* time-domain signal samples.
*
* @param sampleRate
* - sample rate (Hz)
* @param targetFrequency
* - target frequency (Hz)
* @param blockSize
* - number of samples per block of processing
* @param window
* - window to apply against the samples prior to calculating
* the magnitude
*/
public GoertzelFilter( int sampleRate,
int targetFrequency,
int blockSize,
WindowType window )
{
mSampleRate = sampleRate;
mTargetFrequency = targetFrequency;
mBlockSize = blockSize;
mWindowCoefficients = Window.getWindow( window, blockSize );
init();
}
/**
* Establish the pre-calculated values to use in the filter
*/
private void init()
{
double normalizedFrequency = (double) mTargetFrequency / mSampleRate;
mCoefficient = 2.0D * Math.cos( 2 * Math.PI * normalizedFrequency );
}
/**
* Returns the power (dB) of the target frequency within the block of
* samples as normalized against the bin(0) magnitude value
*
* @param samples
* - array of time-domain samples -- array size must be the same
* size as the windowSize parameter passed upon construction
*
* @return - power measurement in dB
* @throws IllegalArgumentException if the sample array size is not equal
* to the defined block size
*/
public int getPower( float[] samples ) throws IllegalArgumentException
{
// Verify size of samples array against block size
if ( samples.length != mBlockSize )
{
throw new IllegalArgumentException(
"Sample array size does not equal Block Size" );
}
// Apply the window against the samples
samples = Window.apply( mWindowCoefficients, samples );
// Process the samples
double s = 0.0D;
double s_prev = 0.0D;
double s_prev2 = 0.0D;
for ( int x = 0; x < samples.length; x++ )
{
s = samples[x] + ( mCoefficient * s_prev ) - s_prev2;
s_prev2 = s_prev;
s_prev = s;
}
//power = s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2 ;
double magnitude = ( s_prev2 * s_prev2 ) + ( s_prev * s_prev) - (mCoefficient * s_prev * s_prev2);
int binZero = getBinZeroPower( samples );
int power = (int) (20 * Math.log10( magnitude / binZero ) );
return power;
}
/**
* Sums all of the sample values to derive the bin 0 power
*/
private int getBinZeroPower( float[] samples )
{
int retVal = 0;
for( int x = 0; x < samples.length; x++ )
{
retVal += samples[ x ];
}
return retVal;
}
/**
* @return the Sample Rate
*/
public int getSampleRate()
{
return mSampleRate;
}
/**
* @return the Target Frequency
*/
public int getTargetFrequency()
{
return mTargetFrequency;
}
/**
* @return the Block Size
*/
public int getBlockSize()
{
return mBlockSize;
}
}

View File

@@ -0,0 +1,197 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
/***************************************************************************
* Copyright (C) 2011 by Paul Lutus *
* lutusp@arachnoid.com *
* *
* 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 2 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
// http://en.wikipedia.org/wiki/Digital_biquad_filter
final public class IIRBiQuadraticFilter {
public enum Type {
BANDPASS, LOWPASS, HIGHPASS, NOTCH, PEAK, LOWSHELF, HIGHSHELF
};
double a0, a1, a2, b0, b1, b2;
double x1, x2, y, y1, y2;
double gain_abs;
Type type;
double center_freq, sample_rate, Q, gainDB;
public IIRBiQuadraticFilter() {
}
public IIRBiQuadraticFilter(Type type, double center_freq, double sample_rate, double Q, double gainDB) {
configure(type, center_freq, sample_rate, Q, gainDB);
}
// constructor without gain setting
public IIRBiQuadraticFilter(Type type, double center_freq, double sample_rate, double Q) {
configure(type, center_freq, sample_rate, Q, 0);
}
public void reset() {
x1 = x2 = y1 = y2 = 0;
}
public double frequency() {
return center_freq;
}
public void configure(Type type, double center_freq, double sample_rate, double Q, double gainDB) {
reset();
Q = (Q == 0) ? 1e-9 : Q;
this.type = type;
this.sample_rate = sample_rate;
this.Q = Q;
this.gainDB = gainDB;
reconfigure(center_freq);
}
public void configure(Type type, double center_freq, double sample_rate, double Q) {
configure(type, center_freq, sample_rate, Q, 0);
}
// allow parameter change while running
public void reconfigure(double cf) {
center_freq = cf;
// only used for peaking and shelving filters
gain_abs = Math.pow(10, gainDB / 40);
double omega = 2 * Math.PI * cf / sample_rate;
double sn = Math.sin(omega);
double cs = Math.cos(omega);
double alpha = sn / (2 * Q);
double beta = Math.sqrt(gain_abs + gain_abs);
switch (type) {
case BANDPASS:
b0 = alpha;
b1 = 0;
b2 = -alpha;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case LOWPASS:
b0 = (1 - cs) / 2;
b1 = 1 - cs;
b2 = (1 - cs) / 2;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case HIGHPASS:
b0 = (1 + cs) / 2;
b1 = -(1 + cs);
b2 = (1 + cs) / 2;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case NOTCH:
b0 = 1;
b1 = -2 * cs;
b2 = 1;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case PEAK:
b0 = 1 + (alpha * gain_abs);
b1 = -2 * cs;
b2 = 1 - (alpha * gain_abs);
a0 = 1 + (alpha / gain_abs);
a1 = -2 * cs;
a2 = 1 - (alpha / gain_abs);
break;
case LOWSHELF:
b0 = gain_abs * ((gain_abs + 1) - (gain_abs - 1) * cs + beta * sn);
b1 = 2 * gain_abs * ((gain_abs - 1) - (gain_abs + 1) * cs);
b2 = gain_abs * ((gain_abs + 1) - (gain_abs - 1) * cs - beta * sn);
a0 = (gain_abs + 1) + (gain_abs - 1) * cs + beta * sn;
a1 = -2 * ((gain_abs - 1) + (gain_abs + 1) * cs);
a2 = (gain_abs + 1) + (gain_abs - 1) * cs - beta * sn;
break;
case HIGHSHELF:
b0 = gain_abs * ((gain_abs + 1) + (gain_abs - 1) * cs + beta * sn);
b1 = -2 * gain_abs * ((gain_abs - 1) + (gain_abs + 1) * cs);
b2 = gain_abs * ((gain_abs + 1) + (gain_abs - 1) * cs - beta * sn);
a0 = (gain_abs + 1) - (gain_abs - 1) * cs + beta * sn;
a1 = 2 * ((gain_abs - 1) - (gain_abs + 1) * cs);
a2 = (gain_abs + 1) - (gain_abs - 1) * cs - beta * sn;
break;
}
// prescale flter constants
b0 /= a0;
b1 /= a0;
b2 /= a0;
a1 /= a0;
a2 /= a0;
}
// provide a static amplitude result for testing
public double result(double f) {
double phi = Math.pow((Math.sin(2.0 * Math.PI * f / (2.0 * sample_rate))), 2.0);
return (Math.pow(b0 + b1 + b2, 2.0) - 4.0 * (b0 * b1 + 4.0 * b0 * b2 + b1 * b2) * phi + 16.0 * b0 * b2 * phi * phi) / (Math.pow(1.0 + a1 + a2, 2.0) - 4.0 * (a1 + 4.0 * a2 + a1 * a2) * phi + 16.0 * a2 * phi * phi);
}
// provide a static decibel result for testing
public double log_result(double f) {
double r;
try {
r = 10 * Math.log10(result(f));
} catch (Exception e) {
r = -100;
}
if(Double.isInfinite(r) || Double.isNaN(r)) {
r = -100;
}
return r;
}
// return the constant set for this filter
public double[] constants() {
return new double[]{b0, b1, b2, a1, a2};
}
// perform one filtering step
public double filter(double x) {
y = b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;
x2 = x1;
x1 = x;
y2 = y1;
y1 = y;
return (y);
}
}

View File

@@ -0,0 +1,51 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
import sample.Listener;
public class SquaringFilter implements Listener<Float>
{
private Listener<Boolean> mListener;
/**
* Transforms a short sample waveform into a square wave
*/
public SquaringFilter()
{
}
@Override
public void receive( Float sample )
{
if( mListener != null )
{
mListener.receive( ( sample >= 0 ? true : false ) );
}
}
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
public void clearListener( Listener<Boolean> listener )
{
mListener = null;
}
}

239
src/dsp/filter/Window.java Normal file
View File

@@ -0,0 +1,239 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.filter;
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
*
*/
public class Window
{
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 );
}
}
public static double[] getRectangularWindow( int length )
{
double[] coefficients = new double[ length ];
Arrays.fill( coefficients, 1.0d );
return coefficients;
}
public static double[] getCosineWindow( int length )
{
double[] coefficients = new double[ length ];
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;
}
public static double[] getBlackmanWindow( 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 ) ) ) );
}
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 ) ) );
}
}
return coefficients;
}
public static double[] getHanningWindow( int 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 ) ) );
}
}
return coefficients;
}
/**
* Apply the window against an array of float-type samples
*/
public static float[] apply( double[] coefficients, float[] samples )
{
for( int x = 0; x < samples.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 < samples.length; x++ )
{
samples[ x ] = samples[ x ] * coefficients[ x ];
}
return 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 ] );
if( breaks )
{
sb.append( "\n" );
}
else
{
if( x % 8 == 7 )
{
sb.append( "\n" );
}
else
{
sb.append( "\t" );
}
}
}
return sb.toString();
}
public static void main( String[] args )
{
System.out.println( arrayToString( getBlackmanWindow( 31 ), false ) );
}
}

View File

@@ -0,0 +1,287 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.fsk;
import instrument.Instrumentable;
import instrument.tap.Tap;
import instrument.tap.stream.SymbolEventTap;
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
import sample.Listener;
import buffer.AveragingCircularBooleanBuffer;
import dsp.fsk.SymbolEvent.Shift;
/**
* Binary Frequency Shift Keying (FSK) decoder. Implements a BFSK correlation
* decoder. Provides normal or inverted decoded output. Requires the symbol
* rate to be an integer multiple of the sample rate.
*
* Implements a correlation decoder that correlates a modulated FSK signal upon
* against a one-baud delayed version of itself, converting all float samples to
* binary (0 or 1), and uses XOR to produce the correlation value. Low pass
* filtering smoothes the correlated output.
*
* Automatically aligns to symbol timing during symbol transitions (0 to 1, or
* 1 to 0) by inspecting the samples at the symbol edges and advancing or
* retarding symbol window to maintain continuous symbol alignment.
*
* Use a DC-removal filter prior to this decoder to ensure samples don't have
* a DC component.
*
* Implements instrumentable interface, so that slice events can be received
* externally to analyze decoder performance.
*/
public class FSK2Decoder implements Instrumentable, Listener<Float>
{
public enum Output{ NORMAL, INVERTED };
/* Instrumentation taps */
private static final String INSTRUMENT_DECISION =
"Tap Point: FSK2 Symbol Decision";
private ArrayList<Tap> mAvailableTaps;
private ArrayList<SymbolEventTap> mTaps = new ArrayList<SymbolEventTap>();
private Listener<Boolean> mListener;
private AveragingCircularBooleanBuffer mDelayBuffer;
private AveragingCircularBooleanBuffer mLowPassFilter;
private Slicer mSlicer;
private boolean mNormalOutput;
private int mSamplesPerSymbol;
private int mSymbolRate;
public FSK2Decoder( int sampleRate, int symbolRate, Output output )
{
/* Ensure we're using an integral symbolRate of the sampleRate */
assert( sampleRate % symbolRate == 0 );
mSamplesPerSymbol = (int)( sampleRate / symbolRate );
mNormalOutput = ( output == Output.NORMAL );
mSymbolRate = symbolRate;
mDelayBuffer = new AveragingCircularBooleanBuffer( mSamplesPerSymbol );
mLowPassFilter = new AveragingCircularBooleanBuffer( mSamplesPerSymbol );
mSlicer = new Slicer( mSamplesPerSymbol );
}
/**
* Disposes of all references to prepare for garbage collection
*/
public void dispose()
{
mListener = null;
}
/**
* Primary sample input
*/
@Override
public void receive( Float floatSample )
{
/* Square the sample. Greater than zero is a 1 (true) and less than
* zero is a 0 (false) */
boolean bitSample = ( floatSample >= 0.0f );
/* Feed the delay buffer and fetch the one-baud delayed sample */
boolean delayedBitSample = mDelayBuffer.get( bitSample );
/* Correlation: xor current bit with delayed bit */
boolean softBit = bitSample ^ delayedBitSample;
/* Low pass filter to smooth the correlated values */
boolean filteredSoftBit = mLowPassFilter.getAverage( softBit );
/* Send the filtered correlated bit to the slicer */
mSlicer.receive( filteredSoftBit );
}
/**
* Registers a listener to receive the decoded FSK bits
*/
public void setListener( Listener<Boolean> listener )
{
mListener = listener;
}
/**
* Removes the listener
*/
public void removeListener( Listener<Boolean> listener )
{
mListener = null;
}
/**
* Symbol slicer with auto-aligning baud timing
*/
public class Slicer
{
private BitSet mBitSet = new BitSet();
private int mSymbolLength;
private int mDecisionThreshold;
private int mSampleCounter;
public Slicer( int samplesPerSymbol )
{
mSymbolLength = samplesPerSymbol;
mDecisionThreshold = (int)( mSymbolLength / 2 );
/* Adjust for an odd number of samples per baud */
if( mSymbolLength % 2 == 1 )
{
mDecisionThreshold++;
}
}
public void receive( boolean softBit )
{
if( mSampleCounter >= 0 )
{
if( softBit )
{
mBitSet.set( mSampleCounter );
}
else
{
mBitSet.clear( mSampleCounter );
}
}
mSampleCounter++;
if( mSampleCounter >= mSymbolLength )
{
boolean decision = mBitSet.cardinality() >= mDecisionThreshold;
send( decision );
/* Shift timing left if the left bit in the bitset is opposite
* the decision and the right bit is the same */
if( ( mBitSet.get( 0 ) ^ decision ) &&
( !( mBitSet.get( mSymbolLength - 1 ) ^ decision ) ) )
{
sendTapEvent( mBitSet, Shift.LEFT, decision );
reset();
mSampleCounter--;
}
/* Shift timing right if the left bit is the same as the
* decision and the right bit is opposite */
else if( ( !( mBitSet.get( 0 ) ^ decision ) ) &&
( mBitSet.get( mSymbolLength - 1 ) ^ decision ) )
{
sendTapEvent( mBitSet, Shift.RIGHT, decision );
/* Last bit from previous symbol to pre-fill next symbol */
boolean previousSoftBit = mBitSet.get( mSymbolLength - 1 );
reset();
if( previousSoftBit )
{
mBitSet.set( 0 );
}
mSampleCounter++;
}
/* No shift */
else
{
sendTapEvent( mBitSet, Shift.NONE, decision );
reset();
}
}
}
/**
* Sends the bit decision to the listener
*/
private void send( boolean decision )
{
if( mListener != null )
{
mListener.receive( mNormalOutput ? decision : !decision );
}
}
private void reset()
{
mBitSet.clear();
mSampleCounter = 0;
}
/**
* Sends instrumentation tap event to all registered listeners
*/
private void sendTapEvent( BitSet bitset, Shift shift, boolean decision )
{
for( SymbolEventTap tap: mTaps )
{
SymbolEvent event =
new SymbolEvent( bitset.get( 0, mSymbolLength ),
mSymbolLength,
decision,
shift );
tap.receive( event );
}
}
}
/**
* Get instrumentation taps
*/
@Override
public List<Tap> getTaps()
{
if( mAvailableTaps == null )
{
mAvailableTaps = new ArrayList<Tap>();
mAvailableTaps.add(
new SymbolEventTap( INSTRUMENT_DECISION, 0, .025f ) );
}
return mAvailableTaps;
}
/**
* Add instrumentation tap
*/
@Override
public void addTap( Tap tap )
{
if( tap instanceof SymbolEventTap && !mTaps.contains( tap ) )
{
mTaps.add( (SymbolEventTap)tap );
}
}
/**
* Remove instrumentation tap
*/
@Override
public void removeTap( Tap tap )
{
mTaps.remove( tap );
}
}

View File

@@ -0,0 +1,75 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.fsk;
import java.util.BitSet;
public class SymbolEvent
{
private BitSet mBitset;
private int mSamplesPerSymbol;
private boolean mDecision;
private Shift mShift;
public SymbolEvent( BitSet bitset, int samplesPerSymbol, boolean decision, Shift shift )
{
mBitset = bitset;
mSamplesPerSymbol = samplesPerSymbol;
mDecision = decision;
mShift = shift;
}
public BitSet getBitSet()
{
return mBitset;
}
public int getSamplesPerSymbol()
{
return mSamplesPerSymbol;
}
public boolean getDecision()
{
return mDecision;
}
public Shift getShift()
{
return mShift;
}
public enum Shift
{
LEFT( "<" ),
RIGHT( ">" ),
NONE( "=" );
private String mLabel;
private Shift( String label )
{
mLabel = label;
}
public String getLabel()
{
return mLabel;
}
}
}

View File

@@ -0,0 +1,100 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.nbfm;
import sample.Listener;
import sample.complex.ComplexSample;
public class FMDiscriminator implements Listener<ComplexSample>
{
private Listener<Float> mListener;
private ComplexSample mPreviousSample = new ComplexSample( 0.0f, 0.0f );
private double mGain;
/**
* Implements a polar discriminator with ArcTangent angle estimator to
* demodulate complex sampled frequency modulated signals.
*
* @param gain - gain to be applied to the demodulated output - can be
* dynamically applied to continuously adjust audio output
*/
public FMDiscriminator( double gain )
{
mGain = gain;
}
public void dispose()
{
mListener = null;
mPreviousSample = null;
}
public synchronized void setGain( double gain )
{
mGain = gain;
}
@Override
public void receive( ComplexSample currentSample )
{
/**
* Multiply the current sample against the complex conjugate of the
* previous sample to derive the phase delta between the two samples
*
* Negating the previous sample quadrature produces the conjugate
*/
double i = ( currentSample.inphase() * mPreviousSample.inphase() ) -
( currentSample.quadrature() * -mPreviousSample.quadrature() );
double q = ( currentSample.quadrature() * mPreviousSample.inphase() ) +
( currentSample.inphase() * -mPreviousSample.quadrature() );
double angle;
//Check for divide by zero
if( i == 0 )
{
angle = 0.0;
}
else
{
/**
* Use the arcus tangent of imaginary (q) divided by real (i) to
* get the phase angle (+/-) which was directly manipulated by the
* original message waveform during the modulation. This value now
* serves as the instantaneous amplitude of the demodulated signal
*/
double denominator = 1.0d / i;
angle = Math.atan( (double)q * denominator );
}
if( mListener != null )
{
mListener.receive( (float)( angle * mGain ) );
}
/**
* Store the current sample to use during the next iteration
*/
mPreviousSample = currentSample;
}
public void setListener( Listener<Float> listener )
{
mListener = listener;
}
}

View File

@@ -0,0 +1,92 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.nbfm;
import sample.Broadcaster;
import sample.Listener;
import sample.complex.ComplexSample;
import dsp.filter.ComplexFIRFilter;
import dsp.filter.FilterFactory;
import dsp.filter.FloatFIRFilter;
import dsp.filter.Window.WindowType;
public class FilteringNBFMDemodulator implements Listener<ComplexSample>
{
private Broadcaster<Float> mBroadcaster = new Broadcaster<Float>();
private ComplexFIRFilter mIQFilter;
private FloatFIRFilter mAudioFilter;
private FMDiscriminator mDiscriminator;
public FilteringNBFMDemodulator( double[] iqFilter, double iqGain,
double[] audioFilter, double audioGain )
{
mIQFilter = new ComplexFIRFilter( iqFilter, iqGain );
mAudioFilter = new FloatFIRFilter( audioFilter, audioGain );
/**
* Add the FM polar discriminator with a gain of 32768, to match
* the fact that we're producing demodulated float values
*/
mDiscriminator = new FMDiscriminator( 32768 );
mIQFilter.setListener( mDiscriminator );
mDiscriminator.setListener( mAudioFilter );
mAudioFilter.setListener( mBroadcaster );
}
/**
* Implements a quadrature narrow-band demodulator that produces float
* valued demodulated, audio filtered output samples.
*/
public FilteringNBFMDemodulator()
{
this( FilterFactory.getLowPass( 48000, 4000, 6000, 48, WindowType.HAMMING, true ),
1.0002,
FilterFactory.getLowPass( 48000, 3000, 6000, 48, WindowType.HAMMING, true ),
2 );
}
/**
* Receive method for complex samples that are fed into this class to be
* processed
*/
@Override
public void receive( ComplexSample quadratureSample )
{
mIQFilter.receive( quadratureSample );
}
/**
* Adds a listener to receive demodulated samples
*/
public void addListener( Listener<Float> listener )
{
mBroadcaster.addListener( listener );
}
/**
* Removes a listener from receiving demodulated samples
*/
public void removeListener( Listener<Float> listener )
{
mBroadcaster.removeListener( listener );
}
}

View File

@@ -0,0 +1,92 @@
/*******************************************************************************
* 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/>
******************************************************************************/
package dsp.nbfm;
import sample.Broadcaster;
import sample.Listener;
import sample.complex.ComplexSample;
import dsp.filter.ComplexFIRFilter;
import dsp.filter.FilterFactory;
import dsp.filter.Window.WindowType;
public class NBFMDemodulator implements Listener<ComplexSample>
{
private Broadcaster<Float> mBroadcaster = new Broadcaster<Float>();
private ComplexFIRFilter mIQFilter;
private FMDiscriminator mDiscriminator;
public NBFMDemodulator( double[] iqFilter, double iqGain )
{
mIQFilter = new ComplexFIRFilter( iqFilter, iqGain );
/**
* Add the FM polar discriminator with a gain of 32768, to match
* the fact that we're producing demodulated float values
*/
mDiscriminator = new FMDiscriminator( 32768 );
mIQFilter.setListener( mDiscriminator );
mDiscriminator.setListener( mBroadcaster );
}
public void dispose()
{
mDiscriminator.dispose();
mBroadcaster.dispose();
mIQFilter.dispose();
}
/**
* Implements a quadrature narrow-band demodulator that produces float
* valued demodulated, audio filtered output samples.
*/
public NBFMDemodulator()
{
this( FilterFactory.getLowPass( 48000, 4000, 73, WindowType.HAMMING ),
1.0002 );
}
/**
* Receive method for complex samples that are fed into this class to be
* processed
*/
@Override
public void receive( ComplexSample sample )
{
mIQFilter.receive( sample );
}
/**
* Adds a listener to receive demodulated samples
*/
public void addListener( Listener<Float> listener )
{
mBroadcaster.addListener( listener );
}
/**
* Removes a listener from receiving demodulated samples
*/
public void removeListener( Listener<Float> listener )
{
mBroadcaster.removeListener( listener );
}
}