/** * Returns the product x * y. * @param x a complex number. * @param y a complex number. * @return the product. */ public static Cdouble mul(Cdouble x, Cdouble y) { return x.times(y); }
/** * Evaluates H(1/z). */ private Cdouble hi( Cdouble z, Cdouble[] poles, Cdouble[] zeros, double gain) { Cdouble c1 = new Cdouble(1.0,0.0); Cdouble hz = new Cdouble(c1); for (int iz=0; iz<_nz; ++iz) { Cdouble zi = zeros[iz]; hz.timesEquals(c1.minus(zi.times(z))); } Cdouble hp = new Cdouble(c1); for (int ip=0; ip<_np; ++ip) { Cdouble pi = poles[ip]; hp.timesEquals(c1.minus(pi.times(z))); } return hz.over(hp).times(gain); }
/** * Evaluates residue of G(z) for the n'th derivative and j'th pole. */ private Cdouble gr(int nd, Cdouble polej, Cdouble[] poles, double gain) { Cdouble pj = polej; Cdouble qj = pj.inv(); Cdouble c1 = new Cdouble(1.0,0.0); Cdouble gz = new Cdouble(c1); if (nd==1) { gz.timesEquals(c1.minus(qj)); gz.timesEquals(c1.plus(pj)); gz.timesEquals(pj); gz.timesEquals(0.5); } else if (nd==2) { gz.timesEquals(c1.minus(qj)); gz.timesEquals(c1.minus(pj)); gz.timesEquals(-1.0); } Cdouble gp = new Cdouble(c1); int np = poles.length; for (int ip=0; ip<np; ++ip) { Cdouble pi = poles[ip]; if (!pi.equals(pj) && !pi.equals(pj.conj())) gp.timesEquals(c1.minus(pi.times(qj))); gp.timesEquals(c1.minus(pi.times(pj))); } return gz.over(gp).times(gain); }
/** * Constructs a recursive 2nd-order filter from poles, zeros, and gain. * The poles must be real or conjugate pairs; likewise for the zeros. * @param pole1 the 1st pole. * @param pole2 the 2nd pole. * @param zero1 the 1st zero. * @param zero2 the 2nd zero. * @param gain the filter gain. */ public Recursive2ndOrderFilter( Cdouble pole1, Cdouble pole2, Cdouble zero1, Cdouble zero2, double gain) { Check.argument(pole1.i==0.0 && pole2.i==0.0 || pole2.r==pole1.r && -pole2.i==pole1.i, "poles are real or conjugate pair"); Check.argument(zero1.i==0.0 && zero2.i==0.0 || zero2.r==zero1.r && -zero2.i==zero1.i, "zeros are real or conjugate pair"); _b0 = (float)(gain); _b1 = (float)(-(zero1.r+zero2.r)*gain); _b2 = (float)((zero1.times(zero2)).r*gain); _a1 = (float)(-(pole1.r+pole2.r)); _a2 = (float)((pole1.times(pole2)).r); }
/** * Evaluates residue of H(z) for the j'th pole. The residue is H(z) * evaluated without division by the factor corresponding to the j'th * pole of H(z), which would be zero. If that pole is complex (has a * non-zero imaginary part), then division by the factor corresponding * to its conjugate mate pole is omitted as well. */ private Cdouble hr( Cdouble polej, Cdouble[] poles, Cdouble[] zeros, double gain) { Cdouble pj = polej; Cdouble qj = pj.inv(); Cdouble c1 = new Cdouble(1.0,0.0); Cdouble hz = new Cdouble(c1); for (int iz=0; iz<_nz; ++iz) { Cdouble zi = zeros[iz]; hz.timesEquals(c1.minus(zi.times(qj))); } Cdouble hp = new Cdouble(c1); for (int ip=0; ip<_np; ++ip) { Cdouble pi = poles[ip]; if (!pi.equals(pj) && !pi.equals(pj.conj())) hp.timesEquals(c1.minus(pi.times(qj))); } return hz.over(hp).times(gain); }
/** * Sets the specified array with a specified linear ramp. * Array values are ca+i1*cb1+i2*cb2. * @param ca value of the first element. * @param cb1 gradient in 1st dimension. * @param cb2 gradient in 2nd dimension. * @param cx the array. */ public static void cramp(Cdouble ca, Cdouble cb1, Cdouble cb2, double[][] cx) { int n2 = cx.length; for (int i2=0; i2<n2; ++i2) cramp(ca.plus(cb2.times((double)i2)),cb1,cx[i2]); }
/** * Sets the specified array with a specified linear ramp. * Array values are ca+i1*cb1+i2*cb2+i3*cb3. * @param ca value of the first element. * @param cb1 gradient in 1st dimension. * @param cb2 gradient in 2nd dimension. * @param cb3 gradient in 3rd dimension. * @param cx the array. */ public static void cramp( Cdouble ca, Cdouble cb1, Cdouble cb2, Cdouble cb3, double[][][] cx) { int n3 = cx.length; for (int i3=0; i3<n3; ++i3) cramp(ca.plus(cb3.times((double)i3)),cb1,cb2,cx[i3]); }
/** * Returns x to the y'th power. * @param x a complex number. * @param y a complex number. * @return x to the y'th power. */ public static Cdouble pow(Cdouble x, Cdouble y) { if (x.r==0.0 && x.i==0.0) return new Cdouble(); return exp(y.times(log(x))); }
private static double computeSigma(double sigma, Cdouble[] poles) { int npole = poles.length; double q = sigma/2.0; Cdouble c1 = new Cdouble(1.0); Cdouble cs = new Cdouble(); for (int ipole=0; ipole<npole; ++ipole) { Cdouble pi = poles[ipole]; double a = pow(pi.abs(),-1.0/q); double t = atan2(pi.i,pi.r)/q; Cdouble b = Cdouble.polar(a,t); Cdouble c = c1.minus(b); Cdouble d = c.times(c); cs.plusEquals(b.times(2.0).over(d)); } return sqrt(cs.r); } }
Cdouble hi = hi(pj,poles,zeros,gain); Cdouble hj = hr(pj,poles,zeros,gain); Cdouble hihj = hi.times(hj); double fb0,fb1,fb2,rb0,rb1,rb2,b0,b1,b2,a1,a2; if (pj.i==0.0) { // if pole is real, ...