/** * 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); }
private static double computeGain(Cdouble[] poles) { int npole = poles.length; Cdouble c1 = new Cdouble(1.0,0.0); Cdouble cg = new Cdouble(c1); for (int ipole=0; ipole<npole; ++ipole) { cg.timesEquals(c1.minus(poles[ipole])); } return cg.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); }
Cdouble cone = new Cdouble(1.0,0.0); Cdouble zinv = Cdouble.polar(1.0,-2.0*DBL_PI*f); Cdouble p = new Cdouble(1.0,0.0); for (Cdouble c: _zeros) p.timesEquals(cone.minus(c.times(zinv))); Cdouble q = new Cdouble(1.0,0.0); for (Cdouble d: _poles) q.timesEquals(cone.minus(d.times(zinv))); Cdouble h = p.over(q); af[jf] = (float)h.abs(); pf[jf] = (float)h.arg();
_poles = new Cdouble[np]; _zeros = new Cdouble[np]; Cdouble c1 = new Cdouble(1.0,0.0); Cdouble c2 = new Cdouble(2.0,0.0); Cdouble zj = (lowpass)?c1.neg():c1; Cdouble gain = new Cdouble(c1); for (int j=0,k=np-1; j<np; ++j,--k) { double theta = ftheta+j*dtheta; Cdouble sj = Cdouble.polar(omegac,theta); _zeros[j] = zj; if (j==k) { _poles[j] = (c2.plus(sj)).over(c2.minus(sj)); _poles[j].i = 0.0; } else if (j<k) { _poles[j] = (c2.plus(sj)).over(c2.minus(sj)); _poles[k] = _poles[j].conj(); gain.timesEquals(sj.over(sj.minus(c2))); } else { gain.timesEquals(c2.over(c2.minus(sj)));
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); } }
/** * Sorts array of poles or zeros. After sorting, any complex conjugate * pairs are first in the array, followed by any real poles or zeros. * Also ensures that any complex poles or zeros have conjugate mates. * @return array of sorted poles or zeros. */ private static Cdouble[] sort(Cdouble[] c) { int n = c.length; Cdouble[] cs = new Cdouble[n]; int ns = 0; for (int i=0; i<n; ++i) { if (!c[i].isReal()) { Cdouble cc = c[i].conj(); int j = 0; while (j<n && !cc.equals(c[j])) ++j; Check.argument(j<n,"complex "+c[i]+" has a conjugate mate"); if (i<j) { cs[ns++] = c[i]; cs[ns++] = c[j]; } } } for (int i=0; i<n; ++i) { if (c[i].isReal()) cs[ns++] = c[i]; } return cs; }
/** * Returns the square root of a complex number. * @param x a complex number. * @return the square root. */ public static Cdouble sqrt(Cdouble x) { if (x.r==0.0) { double t = sqrt(0.5*abs(x.i)); return new Cdouble(t,(x.i<0.0)?-t:t); } else { double t = sqrt(2.0*(abs(x)+abs(x.r))); double u = 0.5*t; return (x.r>0.0) ? new Cdouble(u,x.i/t) : new Cdouble(abs(x.i)/t,(x.i<0.0)?-u:u); } }
/** * Returns the sum z + x, where z is this complex number. * @param x a complex number. * @return z + x. */ public Cdouble plus(Cdouble x) { return new Cdouble(this).plusEquals(x); }
/** * Returns the product z * x, where z is this complex number. * @param x a complex number. * @return z * x. */ public Cdouble times(Cdouble x) { return new Cdouble(this).timesEquals(x); }
/** * Returns the complex conjugate of this complex number. * @return the complex conjugate. */ public Cdouble conj() { return new Cdouble(r,-i); }