Class supporting the discrete Fourier transform
*/ public class dft { private int N; // number of "points" in the DFT private Vector data; // Data on which the DFT is performed, // set up in constructor /** Discreat Fourier Transform constructorThe dft constructor is initialized with a Vector of Java point objects. There are a number of ways to regard the x and y values in the point object. In the context of signal processing, the x values would be time, and the y values would be the magnitude of the signal at that time. If a synthetic function, like sin(x) is sampled, the x component is the argument to the function and the y component is the result.
*/ public dft( Vector d ) { N = 0; data = null; if (d != null) { int len = d.size(); if (len > 0) { // check to see if its a Vector of "point" point p = (point)d.elementAt( 0 ); // the Vector length is > 0 && the type is correct N = len; data = d; } } } // dft constructor /** Return the dft point at mThe DFT algorith is an N2 algorithm. For a DFT point at m, there are O(N) calculations. The basic DFT equation for calculating the DFT point m is shown below. In this equation our data values are stored in the array "x". X[n] is one element of this array. "j" is the sqrt(-1) a.k.a the imaginary number i.
N-1 __ \ / x[n](cos((2Pi*n*m)/N) - j*sin((2Pi*n*m)/N)) -- n= 0
This calculation returns a complex value, with a real part calculated from the cosine part of the equation and the imaginary part calculated from the sine part of the equation.
*/ public complex dftPoint( int m ) { final double twoPi = 2 * Math.PI; complex cx = new complex(0.0f,0.0f); if (m >= 0 && m < N) { double R = 0.0; double I = 0.0; // At m == 0 the DFT reduces to the sum of the data if (m == 0) { point p; for (int n = 0; n < N; n++) { p = (point)data.elementAt( n ); R = R + p.y; } // for } else { double x; double scale; point p; for (int n = 0; n < N; n++) { p = (point)data.elementAt( n ); x = p.y; scale = (twoPi * n * m)/N; R = R + x * Math.cos( scale ); I = I - x * Math.sin( scale ); } // for } // else cx.real( (float)R ); cx.imag( (float)I ); } return cx; } // dftPoint } // class dft dft/util/ 40755 1040 1001 0 7420173035 10707 5 ustar everyone dft/util/complex.java 100644 1040 1001 7257 7420100527 13325 0 ustar everyone package util; import java.lang.Math.*; /** Container for complex numbers */ public class complex { private float R; // real part private float I; // imaginary part /** Argumentless constructor */ public complex() { R = 0.0f; I = 0.0f; } /** constructor with real and imaginary parts */ public complex(float r, float i) { R = r; I = i; } /** constructor with a phasor argument */ public complex( phasor ph ) { R = ph.cx().R; I = ph.cx().I; } /** set the real part */ public void real( float r) { R = r; } /** return the real part */ public float real() { return R; } /** set the imaginary part */ public void imag(float i) { I = i; } /** return the imaginary part */ public float imag() { return I; } /** Return the magnitude of the complex numberThe real and imaginary parts of a complex number are two vectors at right angles to each other. The magnitude of the complex number is the hypotenuse formed by these right angles. So the magnitude is the distance formula:
M = sqrt( R2 + I2 )*/ public float mag() { float Mag = (float)Math.sqrt( (R*R) + (I*I) ); return Mag; } /** Return the angle of the "phasor", in degrees.
The "phasor" is the vector formed from the real and the imaginary parts of the complex number. This vector has a magnitude (returned by mag() above) and an angle returned by this function.
*/ public float angle() { // This result will be in radians float Angle = (float)Math.atan( I / R ); // convert radians to degrees Angle = (Angle * 180.0f)/(float)Math.PI; return Angle; } // angle /** return the "phasor" equivalent to the complex number */ public phasor phase() { phasor ph = new phasor( mag(), angle() ); return ph; } /** Add the rhs complex number to "this" complex number. */ public complex add( complex rhs ) { float newR, newI; newR = R + rhs.R; newI = I + rhs.I; complex C = new complex(newR, newI); return C; } /** Subtract the rhs complex number from "this" complex number. */ public complex sub( complex rhs ) { float newR, newI; newR = R - rhs.R; newI = I - rhs.I; complex C = new complex(newR, newI); return C; } public complex mult( complex rhs ) { float newR, newI; newR = (R * rhs.R) - (I * rhs.I); newI = (R * rhs.I) + (rhs.R * I); complex C = new complex(newR, newI); return C; } public complex mult( float rhs ) { float newR = R * rhs; float newI = rhs * I; complex C = new complex(newR, newI); return C; } /** Divide "this" by the divisor argument */ public complex div( complex divisor ) { float R2 = divisor.R; float I2 = divisor.I; float newR, newI; float div; div = (R2*R2) + (I2*I2); newR = ((R*R2) + (I*I2)) / div; newI = ((R2*I) - (R*I2)) / div; complex C = new complex(newR, newI); return C; } public complex div( float divisor ) { float R2 = divisor; float I2 = 0; float newR, newI; float div; div = R2 * R2; newR = (R * R2) / div; newI = (R2 * I) / div; complex C = new complex(newR, newI); return C; } /** compare rhs to "this". If rhs and this are equal, return true */ public boolean eq( complex rhs ) { return (R == rhs.R && I == rhs.I); } } dft/util/phasor.java 100644 1040 1001 4440 7420100527 13141 0 ustar everyone package util; import java.lang.Math.*; /**A container for magnitude and angle (which complements the complex container)
This class follows the naming convention in Richard Lyon's stellar book Understanding Digital Signal Processing. A "phasor" is sometimes called a vector. However, the term vector also refers to a one dimensional matrix, so the term "phasor". The phasor consists of a magnitude and an angle. It is complementary to a complex number, since a complex number can be formed from a phasor and a phasor can be formed form a complex number.
*/ public class phasor { private float M; // magnitude private float A; // Angle of the magnitude vector, in degrees public phasor() { M = 0.0f; A = 0.0f; } /** Construct a phasor object from a magnitude and an angle. The angle is in degrees */ public phasor( float m, float a ) { M = m; A = a; } /** Construct a phasor object from a complex object. */ public phasor( complex c ) { M = c.mag(); A = c.angle(); } /** set the magnitude */ public void mag( float m ) { M = m; } /** return the magnitude */ public float mag() { return M; } /** set the angle of the phasor */ public void angle( float a ) { A = a; } /** return the angle of the phasor */ public float angle() { return A; } /** return the complex form of the phasorThis function uses the trig. identities:
sin(angle) = I / mag
mag * sin(angle) = I
and
cos(angle) = R / mag
mag * cos(angle) = R
Note that the angle stored in the phasor is in degrees. The Java trig functions take angles in radians.
For a more accurate answer, this function computes the real and imaginary parts of the complex number in double precision and then rounds down to a float.
*/ public complex cx() { double R; double I; double radA = (A * Math.PI)/180; R = Math.cos( radA ) * M; I = Math.sin( radA ) * M; complex c = new complex( (float)R, (float)I ); return c; } } dft/util/point.java 100644 1040 1001 306 7420100527 12753 0 ustar everyone package util; /** Package an point on a graph (e.g., x and y values) */ public class point { public point( float aix, float why ) { x = aix; y = why; } public float x, y; } dft/util/signal.java 100644 1040 1001 12550 7420100527 13143 0 ustar everyone package util; import java.lang.Math; /** This class provides a signal generator for testing Fourier transform code. */ public class signal { /** current x */ private float x; /** number of samples for a signal phase */ private float rate; /** "step" (or increment) between samples */ private float step; /** coeficients scaling the result of the sin function (e.g., coef[i]*sin(x) */ private float coef[]; /** scaling factor for the argument to the sin function (e.g., sin( scale[i]*x ) */ private float scale[]; /** offset factor for the argument to the sin function (e.g., sin(x + offset[i]) */ private float offset[]; /**The constructor is passed:
@param r The number of samples for one cycle of the generated signal. @param c An array of coefficients for the result of the sin function (for example: for c = {1.5f, 3.0f}, 1.5sin(x) + 3.0sin(x)) @param s An array of scaling factors for the argument to the sin function (for example: for s = {0.5f, 0.25f}, sin(x/2) + sin(x/4)) @param o An array of soffset factors for the argument to the sin function (for example: for o = {(float)Math.PI/4.0f}, sin( x + pi/4 ))The signal produced by nextVal() is the sum of a set of sin(x) values. The number of sin(x) values is len = max(c.length, s.length).
. len . --- signal(x) = \ . / ci * sin( si * x + oi) . --- . i=0
Here ci scales the result of the sin(x) function, si scales the argument to sin(x) and oi is an offset to the argument of sin(x).
Example:
float sampleRate = 16;
float coef[] = new float[] {1.0f};
float scale[] = new float[] {2.0f, 1.0f, 0.5f};
float offset[] = new float[] { 0.0f, (float)Math.PI/4.0f }
signal sig = new signal( sampleRate, coef, scale, offset );
This will result in the signal
f(x) = sin( 2x ) + sin( x + pi/4 ) + sin(x/2);
Any or all of the coef, scale and offset arguments may be null. If they are supplied, they do not have to be the same size. The shorter array will be padded with appropriate values by the signal object.
If coef, scale and offset arguments are null the signal will be sin(x).
The signal that would be generated from the coef, scale and offset arguments shown above is shown below.
The loop below is used to generate the points in the graph from the range 0...6Pi.
float range = 3.0f * (2.0f * (float)Math.PI);
point p;
float x;
do {
p = sig.nextVal();
data.addElement( p );
x = p.x;
System.out.println(" " + p.x + " " + p.y );
} while (x <= range);
*/
public signal(float r, float c[], float s[], float o[] )
{
x = (float)0.0;
rate = r;
if (c == null) {
coef = new float[] {1.0f};
}
else {
coef = c;
}
if (s == null) {
scale = new float[] {1.0f};
}
else {
scale = s;
}
if (o == null) {
offset = new float[] {0.0f};
}
else {
offset = o;
}
int maxLen;
maxLen = Math.max( coef.length, scale.length );
maxLen = Math.max( offset.length, maxLen );
if (coef.length < maxLen) {
float newCoef[] = new float[ maxLen ];
int i;
for (i = 0; i < coef.length; i++) {
newCoef[i] = coef[i];
}
for (i = coef.length; i < maxLen; i++) {
newCoef[i] = 1.0f;
}
coef = newCoef;
}
if (scale.length < maxLen) {
float newScale[] = new float[ maxLen ];
int i;
for (i = 0; i < scale.length; i++) {
newScale[i] = scale[i];
}
for (i = scale.length; i < maxLen; i++) {
newScale[i] = 1.0f;
}
scale = newScale;
}
if (offset.length < maxLen) {
float newOffset[] = new float[ maxLen ];
int i;
for (i = 0; i < offset.length; i++) {
newOffset[i] = offset[i];
}
for (i = offset.length; i < maxLen; i++) {
newOffset[i] = 0.0f;
}
offset = newOffset;
}
//
// Adjust the sample rate for the scaling of sin(x)
//
float maxScale = scale[0];
for (int i = 1; i < scale.length; i++) {
maxScale = Math.max( maxScale, scale[i] );
}
rate = rate * maxScale;
step = (float)((2.0 * Math.PI)/(double)rate);
} // signal
/**
Get the next signal value
This function returns the next signal value and then increments the counter by step.
*/ public point nextVal() { point p = new point( x, 0.0f); double scaleX; for (int i = 0; i < coef.length; i++) { scaleX = scale[i] * x; p.y = p.y + coef[i] * (float)Math.sin(scaleX + offset[i]); } x = x + step; return p; } // next_t }