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

Revert "wip: adds filter viewer JavaFX utility application for visualizing filters designed through the filter factory."

Browse Source 
This reverts commit 097884337b.
This commit is contained in:
Denny
2016-04-18 18:44:43 -04:00
parent 097884337b
commit 5cd51c3199
13 changed files with 56 additions and 1100 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/dsp/filter/FilterFactory.java
View File

@@ -1,6 +1,6 @@
/*******************************************************************************
* SDR Trunk
* Copyright (C) 2014-2016,2015 Dennis Sheirer
* Copyright (C) 2014,2015 Dennis Sheirer
*
* Root Raised Cosine filter designer:
* Copyright 2002,2007,2008,2012,2013 Free Software Foundation, Inc.
@@ -35,46 +35,9 @@ import dsp.filter.cic.ComplexPrimeCICDecimate;
public class FilterFactory
{
private final static Logger mLog = LoggerFactory.getLogger( FilterFactory.class );
private final static Logger mLog =
LoggerFactory.getLogger( FilterFactory.class );
public static final double TWO_PI = Math.pow( Math.PI, 2.0 );
public static float[] getNewLowPass( int sampleRate, long passFrequency, long stopFrequency,
int attenuation, WindowType window )
{
assert( passFrequency <= sampleRate / 2 );
int length = getTapCount( sampleRate, passFrequency, stopFrequency, attenuation );
if( length % 2 == 0 )
{
length++;
}
float[] coefficients = new float[ length ];
int half = length / 2;
int K = (int)( length * ( 2.0d * (double)passFrequency / (double) sampleRate ));
double scale = 1.0 / (double)length;
coefficients[ half ] = (float)( K / (double)length );
for( int k = 1; k <= half; k++ )
{
float coefficient = (float)( scale * ( Math.sin( Math.PI * k * K / length ) /
Math.sin( Math.PI * k / length ) ) );
coefficients[ k + half ] = coefficient;
coefficients[ half - k ] = coefficient;
}
Window.apply( WindowType.HANNING, coefficients );
return coefficients;
}
/**
* Generates coefficients for a unity-gain, windowed low-pass filter
* @param sampleRate - hertz
@@ -115,7 +78,6 @@ public class FilterFactory
}
//Apply the window against the coefficients
coefficients = Window.apply( window, coefficients );
//Normalize to unity (1) gain

2
src/dsp/filter/GoertzelFilter.java
View File

@@ -32,7 +32,7 @@ public class GoertzelFilter
private int mSampleRate;
private long mTargetFrequency;
private int mBlockSize;
private float[] mWindowCoefficients;
private double[] mWindowCoefficients;
private double mCoefficient;

87
src/dsp/filter/Window.java
View File

@@ -19,9 +19,6 @@ package dsp.filter;
import java.util.Arrays;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
/**
* Window creator
*
@@ -36,9 +33,7 @@ import org.slf4j.LoggerFactory;
*/
public class Window
{
private final static Logger mLog = LoggerFactory.getLogger( Window.class );
public static float[] getWindow( WindowType type, int length )
public static double[] getWindow( WindowType type, int length )
{
switch( type )
{
@@ -56,18 +51,18 @@ public class Window
}
}
public static float[] getRectangularWindow( int length )
public static double[] getRectangularWindow( int length )
{
float[] coefficients = new float[ length ];
double[] coefficients = new double[ length ];
Arrays.fill( coefficients, 1.0f );
Arrays.fill( coefficients, 1.0d );
return coefficients;
}
public static float[] getCosineWindow( int length )
public static double[] getCosineWindow( int length )
{
float[] coefficients = new float[ length ];
double[] coefficients = new double[ length ];
if( length % 2 == 0 ) //Even length
@@ -76,7 +71,7 @@ public class Window
for( int x = -half; x < length / 2 + 1; x++ )
{
coefficients[ x + half ] = (float)Math.cos(
coefficients[ x + half ] = Math.cos(
( (double)x * Math.PI ) / ( (double)length + 1.0d ) );
}
}
@@ -86,7 +81,7 @@ public class Window
for( int x = -half; x < half + 1; x++ )
{
coefficients[ x + half ] = (float)Math.cos(
coefficients[ x + half ] = Math.cos(
( (double)x * Math.PI ) / ( (double)length + 1.0d ) );
}
}
@@ -94,62 +89,62 @@ public class Window
return coefficients;
}
public static float[] getBlackmanWindow( int length )
public static double[] getBlackmanWindow( int length )
{
float[] coefficients = new float[ length ];
double[] coefficients = new double[ length ];
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = (float)( .426591D -
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 ) ) ) ) );
( .076848D * Math.cos( ( Math.PI * 4.0D * (double)x ) / (double)( length - 1 ) ) ) );
}
return coefficients;
}
public static float[] getHammingWindow( int length )
public static double[] getHammingWindow( int length )
{
float[] coefficients = new float[ length ];
double[] coefficients = new double[ length ];
if( length % 2 == 0 ) //Even length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = (float)( .54D - ( .46D *
Math.cos( ( Math.PI * ( 2.0D * (double)x + 1.0d ) ) / (double)( length - 1 ) ) ) );
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 ] = (float)( .54D - ( .46D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) ) );
coefficients[ x ] = .54D - ( .46D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) );
}
}
return coefficients;
}
public static float[] getHanningWindow( int length )
public static double[] getHanningWindow( int length )
{
float[] coefficients = new float[ length ];
double[] coefficients = new double[ length ];
if( length % 2 == 0 ) //Even length
{
for( int x = 0; x < length; x++ )
{
coefficients[ x ] = (float)( .5D - ( .5D *
Math.cos( ( Math.PI * ( 2.0D * (double)x + 1 ) ) / (double)( length - 1 ) ) ) );
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 ] = (float)( .5D - ( .5D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) ) );
coefficients[ x ] = .5D - ( .5D *
Math.cos( ( 2.0D * Math.PI * (double)x ) / (double)( length - 1 ) ) );
}
}
@@ -160,21 +155,41 @@ public class Window
/**
* Apply the window against an array of float-type samples
*/
public static float[] apply( float[] windowCoefficients, float[] samples )
public static float[] apply( double[] coefficients, float[] samples )
{
for( int x = 0; x < windowCoefficients.length; x++ )
for( int x = 0; x < coefficients.length; x++ )
{
samples[ x ] = samples[ x ] * windowCoefficients[ x ];
samples[ x ] = (float)( samples[ x ] * coefficients[ x ] );
}
return samples;
}
public static float[] apply( WindowType type, float[] samples )
{
double[] coefficients = getWindow( type, samples.length );
return apply( coefficients, samples );
}
/**
* Apply the window against an array of float-type samples
*/
public static double[] apply( double[] coefficients, double[] samples )
{
for( int x = 0; x < coefficients.length; x++ )
{
samples[ x ] = samples[ x ] * coefficients[ x ];
}
return samples;
}
public static float[] apply( WindowType type, float[] samples )
public static double[] apply( WindowType type, double[] samples )
{
float[] windowCoefficients = getWindow( type, samples.length );
double[] coefficients = getWindow( type, samples.length );
return apply( windowCoefficients, samples );
return apply( coefficients, samples );
}
@@ -190,7 +205,7 @@ public class Window
/**
* Utility to log the double arrays
*/
public static String arrayToString( float[] array, boolean breaks )
public static String arrayToString( double[] array, boolean breaks )
{
StringBuilder sb = new StringBuilder();
for( int x = 0; x < array.length; x++ )

574
src/filter/Remez.java
View File

@@ -1,574 +0,0 @@
/**************************************************************************
* Parks-McClellan algorithm for FIR filter design (C version)
*-------------------------------------------------
* Copyright (C) 1995 Jake Janovetz (janovetz@coewl.cen.uiuc.edu)
*
* 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
* 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., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*************************************************************************/
package filter;
public class Remez
{
/* Converted to Java by Iain A Robin, June 1998 */
/* changes by Hanns Holger Rutz, Jan 2010:
* needed to make method 'remez' public
*/
public static final int BANDPASS = 1;
public static final int DIFFERENTIATOR = 2;
public static final int HILBERT = 3;
public static final int NEGATIVE = 0;
public static final int POSITIVE = 1;
public static final int GRIDDENSITY = 16;
public static final int MAXITERATIONS = 40;
double pi = Math.PI;
double pi2 = 2.0*pi;
int numtaps;
public Remez()
{
}
/*******************
* createDenseGrid
*=================
* Creates the dense grid of frequencies from the specified bands.
* Also creates the Desired Frequency Response function (d[]) and
* the Weight function (w[]) on that dense grid
*
*
* INPUT:
* ------
* int r - 1/2 the number of filter coefficients
* int numtaps - Number of taps in the resulting filter
* int numband - Number of bands in user specification
* double[] bands - User-specified band edges [2*numband]
* double[] des - Desired response per band [numband]
* double[] weight - Weight per band [numband]
* int symmetry - Symmetry of filter - used for grid check
*
* OUTPUT:
* -------
* int gridSize - Number of elements in the dense frequency grid
* double[] grid - Frequencies (0 to 0.5) on the dense grid [gridSize]
* double[] d - Desired response on the dense grid [gridSize]
* double[] w - Weight function on the dense grid [gridSize]
*******************/
void createDenseGrid(int r, int numtaps, int numband, double bands[],
double des[], double weight[], int gridSize,
double grid[], double d[], double w[],
int symmetry) {
double lowf, highf;
double delf = 0.5/(GRIDDENSITY*r);
// For differentiator, hilbert,
// symmetry is odd and grid[0] = max(delf, band[0])
if ((symmetry == NEGATIVE) && (delf > bands[0]))
bands[0] = delf;
int j=0;
int k;
for (int band=0; band < numband; band++) {
grid[j] = bands[2*band];
lowf = bands[2*band];
highf = bands[2*band + 1];
k = (int)Math.round((highf - lowf)/delf);
for (int i=0; i<k; i++) {
d[j] = des[band];
w[j] = weight[band];
grid[j] = lowf;
lowf += delf;
j++;
}
grid[j-1] = highf;
}
// Similar to above, if odd symmetry, last grid point can't be .5
// - but, if there are even taps, leave the last grid point at .5
if ((symmetry == NEGATIVE) &&
(grid[gridSize-1] > (0.5 - delf)) && ((numtaps % 2) != 0)) {
grid[gridSize-1] = 0.5-delf;
}
}
/********************
* initialGuess
*==============
* Places Extremal Frequencies evenly throughout the dense grid.
*
*
* INPUT:
* ------
* int r - 1/2 the number of filter coefficients
* int gridSize - Number of elements in the dense frequency grid
*
* OUTPUT:
* -------
* int ext[] - Extremal indexes to dense frequency grid [r+1]
********************/
void initialGuess(int r, int ext[], int gridSize) {
for (int i=0; i<=r; i++) ext[i] = i * (gridSize-1) / r;
}
/***********************
* calcParms
*===========
*
*
* INPUT:
* ------
* int r - 1/2 the number of filter coefficients
* int[] ext - Extremal indexes to dense frequency grid [r+1]
* double[] grid - Frequencies (0 to 0.5) on the dense grid [gridSize]
* double[] d - Desired response on the dense grid [gridSize]
* double[] w - Weight function on the dense grid [gridSize]
*
* OUTPUT:
* -------
* double[] ad - 'b' in Oppenheim & Schafer [r+1]
* double[] x - [r+1]
* double[] y - 'C' in Oppenheim & Schafer [r+1]
***********************/
void calcParms(int r, int ext[], double grid[], double d[], double w[],
double ad[], double x[], double y[]) {
double sign, xi, delta, denom, numer;
// Find x[]
for (int i=0; i<=r; i++) x[i] = Math.cos(pi2 * grid[ext[i]]);
// Calculate ad[] - Oppenheim & Schafer eq 7.132
int ld = (r-1)/15 + 1; // Skips around to avoid round errors
for (int i=0; i<=r; i++) {
denom = 1.0;
xi = x[i];
for (int j=0; j<ld; j++) {
for (int k=j; k<=r; k+=ld)
if (k != i) denom *= 2.0*(xi - x[k]);
}
if (Math.abs(denom)<0.00001) denom = 0.00001;
ad[i] = 1.0/denom;
}
// Calculate delta - Oppenheim & Schafer eq 7.131
numer = denom = 0;
sign = 1;
for (int i=0; i<=r; i++) {
numer += ad[i] * d[ext[i]];
denom += sign * ad[i]/w[ext[i]];
sign = -sign;
}
delta = numer/denom;
sign = 1;
// Calculate y[] - Oppenheim & Schafer eq 7.133b
for (int i=0; i<=r; i++) {
y[i] = d[ext[i]] - sign * delta/w[ext[i]];
sign = -sign;
}
}
/*********************
* computeA
*==========
* Using values calculated in CalcParms, ComputeA calculates the
* actual filter response at a given frequency (freq). Uses
* eq 7.133a from Oppenheim & Schafer.
*
*
* INPUT:
* ------
* double freq - Frequency (0 to 0.5) at which to calculate A
* int r - 1/2 the number of filter coefficients
* double[] ad - 'b' in Oppenheim & Schafer [r+1]
* double[] x - [r+1]
* double[] y - 'C' in Oppenheim & Schafer [r+1]
*
* OUTPUT:
* -------
* Returns double value of A[freq]
*********************/
double computeA(double freq, int r, double ad[], double x[], double y[]) {
double c;
double numer = 0;
double denom = 0;
double xc = Math.cos(pi2 * freq);
for (int i=0; i<=r; i++) {
c = xc - x[i];
if (Math.abs(c) < 1.0e-7) {
numer = y[i];
denom = 1;
break;
}
c = ad[i]/c;
denom += c;
numer += c*y[i];
}
return numer/denom;
}
/************************
* calcError
*===========
* Calculates the Error function from the desired frequency response
* on the dense grid (d[]), the weight function on the dense grid (w[]),
* and the present response calculation (A[])
*
*
* INPUT:
* ------
* int r - 1/2 the number of filter coefficients
* double[] ad - [r+1]
* double[] x - [r+1]
* double[] y - [r+1]
* int gridSize - Number of elements in the dense frequency grid
* double[] grid - Frequencies on the dense grid [gridSize]
* double[] d - Desired response on the dense grid [gridSize]
* double[] w - Weight function on the desnse grid [gridSize]
*
* OUTPUT:
* -------
* double[] e - Error function on dense grid [gridSize]
************************/
void calcError(int r, double ad[], double x[], double y[],
int gridSize, double grid[],
double d[], double w[], double e[]) {
double A;
for (int i=0; i<gridSize; i++) {
A = computeA(grid[i], r, ad, x, y);
e[i] = w[i] * (d[i] - A);
}
}
/************************
* search
*========
* Searches for the maxima/minima of the error curve. If more than
* r+1 extrema are found, it uses the following heuristic (thanks
* Chris Hanson):
* 1) Adjacent non-alternating extrema deleted first.
* 2) If there are more than one excess extrema, delete the
* one with the smallest error. This will create a non-alternation
* condition that is fixed by 1).
* 3) If there is exactly one excess extremum, delete the smaller
* of the first/last extremum
*
*
* INPUT:
* ------
* int r - 1/2 the number of filter coefficients
* int[] ext - Indexes to grid[] of extremal frequencies [r+1]
* int gridSize - Number of elements in the dense frequency grid
* double[] e - Array of error values. [gridSize]
*
* OUTPUT:
* -------
* int[] ext - New indexes to extremal frequencies [r+1]
************************/
void search(int r, int ext[], int gridSize, double e[]) {
boolean up, alt;
int[] foundExt = new int[gridSize]; /* Array of found extremals */
int k = 0;
// Check for extremum at 0.
if (((e[0]>0.0) && (e[0]>e[1])) || ((e[0]<0.0) && (e[0]<e[1])))
foundExt[k++] = 0;
// Check for extrema inside dense grid
for (int i=1; i<gridSize-1; i++) {
if (((e[i]>=e[i-1]) && (e[i]>e[i+1]) && (e[i]>0.0)) ||
((e[i]<=e[i-1]) && (e[i]<e[i+1]) && (e[i]<0.0)))
foundExt[k++] = i;
}
// Check for extremum at 0.5
int j = gridSize-1;
if (((e[j]>0.0) && (e[j]>e[j-1])) || ((e[j]<0.0) && (e[j]<e[j-1])))
foundExt[k++] = j;
// Remove extra extremals
int extra = k - (r+1);
int l;
while (extra > 0) {
up = e[foundExt[0]] > 0.0;
// up = true --> first one is a maximum
// up = false --> first one is a minimum
l=0;
alt = true;
for (j=1; j<k; j++) {
if (Math.abs(e[foundExt[j]]) < Math.abs(e[foundExt[l]]))
l = j; // new smallest error.
if (up && (e[foundExt[j]] < 0.0))
up = false; // switch to a minima
else if (!up && (e[foundExt[j]] > 0.0))
up = true; // switch to a maxima
else {
alt = false;
break; // Ooops, found two non-alternating
} // extrema. Delete smallest of them
} // if the loop finishes, all extrema are alternating
// If there's only one extremal and all are alternating,
// delete the smallest of the first/last extremals.
if (alt && (extra == 1)) {
if (Math.abs(e[foundExt[k-1]]) < Math.abs(e[foundExt[0]]))
l = foundExt[k-1]; // Delete last extremal
else
l = foundExt[0]; // Delete first extremal
}
for (j=l; j<k; j++) // Loop that does the deletion
foundExt[j] = foundExt[j+1];
k--;
extra--;
}
// Copy found extremals to ext[]
for (int i=0; i<=r; i++) ext[i] = foundExt[i];
}
/*********************
* freqSample
*============
* Simple frequency sampling algorithm to determine the impulse
* response h[] from A's found in ComputeA
*
*
* INPUT:
* ------
* int N - Number of filter coefficients
* double[] A - Sample points of desired response [N/2]
* int symmetry - Symmetry of desired filter
*
* OUTPUT:
* -------
* double[] h - Impulse Response of final filter [N]
*********************/
double[] freqSample(double A[], int symm) {
double x, val;
int N = numtaps;
double M = (N-1.0)/2.0;
double[] h = new double[N];
if (symm == POSITIVE) {
if (N%2 != 0) {
for (int n=0; n<N; n++) {
val = A[0];
x = pi2 * (n - M)/N;
for (int k=1; k<=M; k++) val += 2.0 * A[k] * Math.cos(x*k);
h[n] = val/N;
}
}
else {
for (int n=0; n<N; n++) {
val = A[0];
x = pi2 * (n - M)/N;
for (int k=1; k<=(N/2-1); k++) val += 2.0 * A[k] * Math.cos(x*k);
h[n] = val/N;
}
}
}
else {
if (N%2 != 0) {
for (int n=0; n<N; n++) {
val = 0;
x = pi2 * (n - M)/N;
for (int k=1; k<=M; k++) val += 2.0 * A[k] * Math.sin(x*k);
h[n] = val/N;
}
}
else {
for (int n=0; n<N; n++) {
val = A[N/2] * Math.sin(pi * (n - M));
x = pi2 * (n - M)/N;
for (int k=1; k<=(N/2-1); k++) val += 2.0 * A[k] * Math.sin(x*k);
h[n] = val/N;
}
}
}
return h;
}
/*******************
* isDone
*========
* Checks to see if the error function is small enough to consider
* the result to have converged.
*
* INPUT:
* ------
* int r - 1/2 the number of filter coeffiecients
* int[] ext - Indexes to extremal frequencies [r+1]
* double[] e - Error function on the dense grid [gridSize]
*
* OUTPUT:
* -------
* Returns true if the result converged
* Returns false if the result has not converged
********************/
boolean isDone(int r, int ext[], double e[]) {
double min, max, current;
min = max = Math.abs(e[ext[0]]);
for (int i=1; i<=r; i++){
current = Math.abs(e[ext[i]]);
if (current < min) min = current;
if (current > max) max = current;
}
if (((max-min)/max) < 0.0001) return true;
return false;
}
/********************
* remez
*=======
* Calculates the optimal (in the Chebyshev/minimax sense)
* FIR filter impulse response given a set of band edges,
* the desired reponse on those bands, and the weight given to
* the error in those bands.
*
* INPUT:
* ------
* int numtaps - Number of filter coefficients
* int numband - Number of bands in filter specification
* double[] bands - User-specified band edges [2 * numband]
* double[] des - User-specified band responses [numband]
* double[] weight - User-specified error weights [numband]
* int type - Type of filter
*
* OUTPUT:
* -------
* double[] h - Impulse response of final filter [numtaps]
********************/
public double[] remez(int n,
double bands[], double des[], double weight[], int type) {
double c;
numtaps = n;
int symmetry = (type == BANDPASS)? POSITIVE : NEGATIVE;
int r = numtaps/2; // number of extrema
if ((numtaps%2 != 0) && (symmetry == POSITIVE)) r++;
// Predict dense grid size in advance for array sizes
int gridSize = 0;
int numband = des.length;
for (int i=0; i<numband; i++) {
gridSize += (int)Math.round(2*r*GRIDDENSITY*(bands[2*i+1] - bands[2*i]));
}
if (symmetry == NEGATIVE) gridSize--;
double[] grid = new double[gridSize];
double[] d = new double[gridSize];
double[] w = new double[gridSize];
double[] e = new double[gridSize];
double[] x = new double[r+1];
double[] y = new double[r+1];
double[] ad = new double[r+1];
double[] taps = new double[r+1];
int[] ext = new int[r+1];
// Create dense frequency grid
createDenseGrid(r, numtaps, numband, bands, des, weight,
gridSize, grid, d, w, symmetry);
initialGuess(r, ext, gridSize);
// For Differentiator: (fix grid)
if (type == DIFFERENTIATOR) {
for (int i=0; i<gridSize; i++) {
if (d[i] > 0.0001) w[i] /= grid[i];
}
}
// For odd or Negative symmetry filters, alter the
// d[] and w[] according to Parks McClellan
if (symmetry == POSITIVE) {
if (numtaps % 2 == 0) {
for (int i=0; i<gridSize; i++) {
c = Math.cos(pi * grid[i]);
d[i] /= c;
w[i] *= c;
}
}
}
else {
if (numtaps % 2 != 0) {
for (int i=0; i<gridSize; i++) {
c = Math.sin(pi2 * grid[i]);
d[i] /= c;
w[i] *= c;
}
}
else {
for (int i=0; i<gridSize; i++) {
c = Math.sin(pi * grid[i]);
d[i] /= c;
w[i] *= c;
}
}
}
// Perform the Remez Exchange algorithm
int iter;
for (iter=0; iter<MAXITERATIONS; iter++) {
calcParms(r, ext, grid, d, w, ad, x, y);
calcError(r, ad, x, y, gridSize, grid, d, w, e);
search(r, ext, gridSize, e);
if (isDone(r, ext, e)) break;
}
if (iter == MAXITERATIONS) {
System.out.println("Reached maximum iteration count.\nResults may be bad.\n");
}
calcParms(r, ext, grid, d, w, ad, x, y);
// Find the 'taps' of the filter for use with Frequency
// Sampling. If odd or Negative symmetry, fix the taps
// according to Parks McClellan
for (int i=0; i<=numtaps/2; i++) {
if (symmetry == POSITIVE) {
if (numtaps%2 != 0) c = 1;
else c = Math.cos(pi * (double)i/numtaps);
}
else {
if (numtaps%2 != 0) c = Math.sin(pi2 * (double)i/numtaps);
else c = Math.sin(pi * (double)i/numtaps);
}
taps[i] = computeA((double)i/numtaps, r, ad, x, y)*c;
}
// Frequency sampling design with calculated taps
return freqSample(taps, symmetry);
}
}

207
src/filter/viewer/CurveFittedAreaChart.java
View File

@@ -1,207 +0,0 @@
package filter.viewer;
import javafx.collections.ObservableList;
import javafx.geometry.Point2D;
import javafx.scene.Group;
import javafx.scene.chart.AreaChart;
import javafx.scene.chart.NumberAxis;
import javafx.scene.chart.XYChart;
import javafx.scene.shape.ClosePath;
import javafx.scene.shape.CubicCurveTo;
import javafx.scene.shape.LineTo;
import javafx.scene.shape.MoveTo;
import javafx.scene.shape.Path;
import javafx.scene.shape.PathElement;
import javafx.util.Pair;
public class CurveFittedAreaChart extends AreaChart<Number, Number>
{
public CurveFittedAreaChart( NumberAxis xAxis, NumberAxis yAxis )
{
super( xAxis, yAxis );
}
@Override
protected void layoutPlotChildren()
{
super.layoutPlotChildren();
for ( int seriesIndex = 0; seriesIndex < getDataSize(); seriesIndex++ )
{
final XYChart.Series<Number, Number> series = getData()
.get( seriesIndex );
final Path seriesLine = (Path) ( (Group) series.getNode() )
.getChildren().get( 1 );
final Path fillPath = (Path) ( (Group) series.getNode() )
.getChildren().get( 0 );
smooth( seriesLine.getElements(), fillPath.getElements() );
}
}
private int getDataSize()
{
final ObservableList<XYChart.Series<Number, Number>> data = getData();
return ( data != null ) ? data.size() : 0;
}
private static void smooth( ObservableList<PathElement> strokeElements,
ObservableList<PathElement> fillElements )
{
// as we do not have direct access to the data, first recreate the list
// of all the data points we have
final Point2D[] dataPoints = new Point2D[ strokeElements.size() ];
for ( int i = 0; i < strokeElements.size(); i++ )
{
final PathElement element = strokeElements.get( i );
if ( element instanceof MoveTo )
{
final MoveTo move = (MoveTo) element;
dataPoints[ i ] = new Point2D( move.getX(), move.getY() );
} else if ( element instanceof LineTo )
{
final LineTo line = (LineTo) element;
final double x = line.getX(), y = line.getY();
dataPoints[ i ] = new Point2D( x, y );
}
}
// next we need to know the zero Y value
final double zeroY = ( (MoveTo) fillElements.get( 0 ) ).getY();
// now clear and rebuild elements
strokeElements.clear();
fillElements.clear();
Pair<Point2D[], Point2D[]> result = calcCurveControlPoints(
dataPoints );
Point2D[] firstControlPoints = result.getKey();
Point2D[] secondControlPoints = result.getValue();
// start both paths
strokeElements.add(
new MoveTo( dataPoints[ 0 ].getX(), dataPoints[ 0 ].getY() ) );
fillElements.add( new MoveTo( dataPoints[ 0 ].getX(), zeroY ) );
fillElements.add(
new LineTo( dataPoints[ 0 ].getX(), dataPoints[ 0 ].getY() ) );
// add curves
for ( int i = 1; i < dataPoints.length; i++ )
{
final int ci = i - 1;
strokeElements
.add( new CubicCurveTo( firstControlPoints[ ci ].getX(),
firstControlPoints[ ci ].getY(),
secondControlPoints[ ci ].getX(),
secondControlPoints[ ci ].getY(),
dataPoints[ i ].getX(), dataPoints[ i ].getY() ) );
fillElements.add( new CubicCurveTo( firstControlPoints[ ci ].getX(),
firstControlPoints[ ci ].getY(),
secondControlPoints[ ci ].getX(),
secondControlPoints[ ci ].getY(), dataPoints[ i ].getX(),
dataPoints[ i ].getY() ) );
}
// end the paths
fillElements.add( new LineTo(
dataPoints[ dataPoints.length - 1 ].getX(), zeroY ) );
fillElements.add( new ClosePath() );
}
/**
* Calculate open-ended Bezier Spline Control Points.
*/
public static Pair<Point2D[], Point2D[]> calcCurveControlPoints(
Point2D[] dataPoints )
{
Point2D[] firstControlPoints;
Point2D[] secondControlPoints;
int n = dataPoints.length - 1;
if ( n == 1 )
{ // Special case: Bezier curve should be a straight line.
firstControlPoints = new Point2D[ 1 ];
// 3P1 = 2P0 + P3
firstControlPoints[ 0 ] = new Point2D(
( 2 * dataPoints[ 0 ].getX() + dataPoints[ 1 ].getX() ) / 3,
( 2 * dataPoints[ 0 ].getY() + dataPoints[ 1 ].getY() )
/ 3 );
secondControlPoints = new Point2D[ 1 ];
// P2 = 2P1 – P0
secondControlPoints[ 0 ] = new Point2D(
2 * firstControlPoints[ 0 ].getX() - dataPoints[ 0 ].getX(),
2 * firstControlPoints[ 0 ].getY()
- dataPoints[ 0 ].getY() );
return new Pair<Point2D[], Point2D[]>( firstControlPoints,
secondControlPoints );
}
// Calculate first Bezier control points
// Right hand side vector
double[] rhs = new double[ n ];
// Set right hand side X values
for ( int i = 1; i < n - 1; ++i )
{
rhs[ i ] = 4 * dataPoints[ i ].getX()
+ 2 * dataPoints[ i + 1 ].getX();
}
rhs[ 0 ] = dataPoints[ 0 ].getX() + 2 * dataPoints[ 1 ].getX();
rhs[ n - 1 ] = ( 8 * dataPoints[ n - 1 ].getX()
+ dataPoints[ n ].getX() ) / 2.0;
// Get first control points X-values
double[] x = GetFirstControlPoints( rhs );
// Set right hand side Y values
for ( int i = 1; i < n - 1; ++i )
{
rhs[ i ] = 4 * dataPoints[ i ].getY()
+ 2 * dataPoints[ i + 1 ].getY();
}
rhs[ 0 ] = dataPoints[ 0 ].getY() + 2 * dataPoints[ 1 ].getY();
rhs[ n - 1 ] = ( 8 * dataPoints[ n - 1 ].getY()
+ dataPoints[ n ].getY() ) / 2.0;
// Get first control points Y-values
double[] y = GetFirstControlPoints( rhs );
// Fill output arrays.
firstControlPoints = new Point2D[ n ];
secondControlPoints = new Point2D[ n ];
for ( int i = 0; i < n; ++i )
{
// First control point
firstControlPoints[ i ] = new Point2D( x[ i ], y[ i ] );
// Second control point
if ( i < n - 1 )
{
secondControlPoints[ i ] = new Point2D(
2 * dataPoints[ i + 1 ].getX() - x[ i + 1 ],
2 * dataPoints[ i + 1 ].getY() - y[ i + 1 ] );
} else
{
secondControlPoints[ i ] = new Point2D(
( dataPoints[ n ].getX() + x[ n - 1 ] ) / 2,
( dataPoints[ n ].getY() + y[ n - 1 ] ) / 2 );
}
}
return new Pair<Point2D[], Point2D[]>( firstControlPoints,
secondControlPoints );
}
/**
* Solves a tridiagonal system for one of coordinates (x or y) of first
* Bezier control points.
*/
private static double[] GetFirstControlPoints( double[] rhs )
{
int n = rhs.length;
double[] x = new double[ n ]; // Solution vector.
double[] tmp = new double[ n ]; // Temp workspace.
double b = 2.0;
x[ 0 ] = rhs[ 0 ] / b;
for ( int i = 1; i < n; i++ )
{// Decomposition and forward substitution.
tmp[ i ] = 1 / b;
b = ( i < n - 1 ? 4.0 : 3.5 ) - tmp[ i ];
x[ i ] = ( rhs[ i ] - x[ i - 1 ] ) / b;
}
for ( int i = 1; i < n; i++ )
{
x[ n - i - 1 ] -= tmp[ n - i ] * x[ n - i ]; // Backsubstitution.
}
return x;
}
}

39
src/filter/viewer/FilterView.css
View File

@@ -1,39 +0,0 @@
/* ....Show License.... */
.chart {
-fx-background-image: url("../../../images/Curve-fitted-background.png");
-fx-padding: 15 25 15 15;
}
.chart-plot-background {
-fx-background-color: red;
-fx-background-image: url("../../../images/Curve-fitted-chart-background.png"),
url("../../../images/Curve-fitted-graph-gridlines.png");
-fx-background-size: cover, auto;
-fx-background-repeat: no-repeat, repeat;
-fx-background-position: 0% 0%, 0% 100%;
-fx-border-color: black black transparent transparent;
}
.chart-area-symbol {
-fx-background-color: white;
}
.chart-series-area-line {
-fx-stroke: white;
-fx-stroke-width: 2px;
}
.chart-series-area-fill {
-fx-fill: linear-gradient(to right, white, rgba(255,255,255,0));
-fx-blend-mode: OVERLAY;
}
.axis {
-fx-tick-mark-visible: false;
-fx-minor-tick-visible: false;
-fx-tick-length: 3;
-fx-minor-tick-length: 0;
-fx-border-color: transparent;
-fx-tick-label-fill: white;
}
.axis:bottom {
-fx-border-color: black transparent transparent transparent;
}
.axis:left {
-fx-border-color: transparent black transparent transparent;
}

168
src/filter/viewer/FilterView.java
View File

@@ -1,168 +0,0 @@
package filter.viewer;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import dsp.filter.FilterFactory;
import dsp.filter.Window.WindowType;
import javafx.collections.ObservableList;
import javafx.scene.chart.NumberAxis;
import javafx.scene.chart.XYChart;
import javafx.scene.chart.XYChart.Series;
import javafx.scene.layout.BorderPane;
import javafx.scene.layout.HBox;
import javafx.scene.layout.Priority;
public class FilterView extends BorderPane
{
private final static Logger mLog = LoggerFactory.getLogger( FilterView.class );
private CurveFittedAreaChart mImpulseChart;
private NumberAxis mImpulseXAxis;
private NumberAxis mImpulseYAxis;
private CurveFittedAreaChart mDFTChart;
private NumberAxis mDFTXAxis;
private NumberAxis mDFTYAxis;
private static final int FFT_SIZE = 4096;
private FloatFFT_1D mFFT = new FloatFFT_1D( FFT_SIZE );
public FilterView()
{
init();
float[] taps = FilterFactory.getNewLowPass( 8000, 1500, 2500, 50, WindowType.HANNING );
setImpulseTaps( taps );
float[] dft = new float[ FFT_SIZE ];
System.arraycopy( taps, 0, dft, 0, taps.length );
mFFT.realForward( dft );
setDFTTaps( convertDFTToDecibel( dft, taps.length ) );
}
private void logTaps( float[] taps )
{
StringBuilder sb = new StringBuilder();
sb.append( "\n" );
for( float tap: taps )
{
sb.append( "Tap: " + tap + "\n" );
}
mLog.debug( sb.toString() );
}
private void init()
{
HBox hbox = new HBox();
mImpulseXAxis = new NumberAxis( 0.0, 0.5, 1.0 );
mImpulseYAxis = new NumberAxis( -0.2, 1.2, .1 );
mImpulseChart = new CurveFittedAreaChart( mImpulseXAxis, mImpulseYAxis );
mImpulseChart.setLegendVisible( false );
mImpulseChart.setHorizontalGridLinesVisible( false );
mImpulseChart.setVerticalGridLinesVisible(false);
mImpulseChart.setAlternativeColumnFillVisible(false);
mImpulseChart.setAlternativeRowFillVisible(false);
String css = FilterView.class.getResource( "FilterView.css" ).toExternalForm();
mImpulseChart.getStylesheets().add( css );
hbox.getChildren().add( mImpulseChart );
HBox.setHgrow( mImpulseChart, Priority.ALWAYS );
mDFTXAxis = new NumberAxis( 0.0, FFT_SIZE / 2, FFT_SIZE / 8 );
mDFTYAxis = new NumberAxis( -200.0, 10.0, 10.0 );
mDFTChart = new CurveFittedAreaChart( mDFTXAxis, mDFTYAxis );
mDFTChart.setLegendVisible( false );
mDFTChart.setHorizontalGridLinesVisible( false );
mDFTChart.setVerticalGridLinesVisible(false);
mDFTChart.setAlternativeColumnFillVisible(false);
mDFTChart.setAlternativeRowFillVisible(false);
mDFTChart.getStylesheets().add( css );
hbox.getChildren().add( mDFTChart );
HBox.setHgrow( mDFTChart, Priority.ALWAYS );
setCenter( hbox );
}
public void setImpulseTaps( float[] taps )
{
double smallest = 0.0;
double largest = 0.0;
final XYChart.Series<Number, Number> series = new XYChart.Series<>();
for( int x = 0; x < taps.length; x++ )
{
if( taps[ x ] < smallest )
{
smallest = taps[ x ];
}
if( taps[ x ] > largest )
{
largest = taps[ x ];
}
series.getData().add( new XYChart.Data<Number, Number>( x, taps[ x ] ) );
}
mImpulseXAxis.setLowerBound( 0.0 );
mImpulseXAxis.setUpperBound( taps.length - 1 );
mImpulseXAxis.setTickUnit( 1.0 );
double padding = Math.abs( largest - smallest ) * .1;
mImpulseYAxis.setLowerBound( smallest - padding );
mImpulseYAxis.setUpperBound( largest + padding );
mImpulseYAxis.setTickUnit( padding );
ObservableList<Series<Number,Number>> list = mImpulseChart.getData();
list.add( series );
}
public void setDFTTaps( float[] taps )
{
double smallest = 0.0;
final XYChart.Series<Number, Number> series = new XYChart.Series<>();
for( int x = 0; x < taps.length; x++ )
{
if( taps[ x ] < smallest )
{
smallest = taps[ x ];
}
series.getData().add( new XYChart.Data<Number, Number>( x, taps[ x ] ) );
}
ObservableList<Series<Number,Number>> list = mDFTChart.getData();
list.add( series );
}
public float[] convertDFTToDecibel( float[] dft, int tapLength )
{
float[] decibels = new float[ dft.length / 2 ];
int index = 0;
for( int x = 0; x < decibels.length; x ++ )
{
index = x * 2;
decibels[ x ] = 20.0f * (float)Math.log10( ( ( dft[ index ] * dft[ index ] ) +
( dft[ index + 1 ] * dft[ index + 1 ] ) ) );
}
return decibels;
}
}

27
src/filter/viewer/FilterViewer.java
View File

@@ -1,27 +0,0 @@
package filter.viewer;
import javafx.application.Application;
import javafx.scene.Scene;
import javafx.stage.Stage;
public class FilterViewer extends Application
{
public FilterViewer()
{
}
@Override
public void start( Stage primaryStage ) throws Exception
{
Scene scene = new Scene( new FilterView() );
primaryStage.setTitle( "Filter Viewer" );
primaryStage.setScene( scene );
primaryStage.show();
}
public static void main( String[] args )
{
launch( args );
}
}

2
src/spectrum/DFTProcessor.java
View File

@@ -72,7 +72,7 @@ public class DFTProcessor implements Listener<ComplexBuffer>,
private DFTSize mDFTSize = DFTSize.FFT04096;
private DFTSize mNewDFTSize = DFTSize.FFT04096;
private float[] mWindow;
private double[] mWindow;
/* The Cosine and Hanning windows seem to offer the best spectral display
* with minimal bin leakage/smearing */