kaiserfilter.java

用Java 写的Kaiser Filter设计

class KaiserFilter {

	public static final int LP = 1;
	public static final int HP = 2;
	public static final int BP = 3;

	int order, filterType, freqPoints;
	float alpha, f1, f2, fN, atten, trband;
	float[] a;

    public void KaiserFilter() {
        // initial settings
        filterType = LP;
        trband = 500.0f;
        atten = 60.0f;
        // following settings based on default 8000 samples/s rate:
        fN = 4000.0f;
        f1 = 0.0f;
        f2 = 1000.0f;
    }

	float sqr (float x) {
	    return x*x;
	}

	float I0 (float x) {
       	// zero order Bessel function of the first kind
		float eps = 1.0e-6f;   // accuracy parameter
		float fact = 1.0f;
		float x2 = 0.5f * x;
		float p = x2;
		float t = p * p;
		float s = 1.0f + t;
		for (int k = 2; t > eps; k++) {
			p *= x2;
			fact *= k;
			t = sqr(p / fact);
			s += t;
		}
		return s;
	}

    public void setFilterType(String ft) {
        if (ft.equals("LP")) filterType = LP;
        if (ft.equals("BP")) filterType = BP;
        if (ft.equals("HP")) filterType = HP;
    }

    public void setFilterType(int ft) {
        filterType = ft;
    }

    public void setRate(float rate) {
        fN = 0.5f * rate;
    }

    public void setFreq1(float f1) {
        this.f1 = f1;
    }

    public void setFreq2(float f2) {
        this.f2 = f2;
    }

    public void setFreqPoints(int fp) {
        freqPoints = fp;
    }

    public void setAtten(float a) {
	    atten = a;
    }

    public void setTrBand(float tb) {
	    trband = tb;
    }

	public float getCoeff(int i) {
	    return a[i];
	}

	public int estimatedOrder() {

	    // estimate filter order
        order = 2 * (int) ((atten - 7.95) / (14.36*trband/fN) + 1.0f);
        // estimate Kaiser window parameter (alpha):
        if (atten >= 50.0f) alpha = 0.1102f*(atten - 8.7f);
		else
	    if (atten > 21.0f)
            alpha = 0.5842f*(float)Math.exp(0.4f*(float)Math.log(atten - 21.0f))
                        + 0.07886f*(atten - 21.0f);
        if (atten <= 21.0f) alpha = 0.0f;
	    return order;
	}

	public void design() {

		// window function values
	    float I0alpha = I0(alpha);
        int m = order / 2;
	    float[] win = new float[m+1];
	    for (int n=1; n <= m; n++)
            win[n] = I0(alpha*(float)Math.sqrt(1.0f - sqr((float)n/m))) / I0alpha;

        float w0 = 0.0f;
        float w1 = 0.0f;
		switch (filterType) {
			case LP: w0 = 0.0f;
			         w1 = (float)Math.PI*(f2 + 0.5f*trband)/fN;
			break;
			case HP: w0 = (float)Math.PI;
				     w1 = (float)Math.PI*(1.0f - (f1 - 0.5f*trband)/fN);
			break;
			case BP: w0 = 0.5f * (float)Math.PI * (f1 + f2) / fN;
				     w1 = 0.5f * (float)Math.PI * (f2 - f1 + trband) / fN;
			break;
		}

	    // filter coefficients (NB not normalised to unit maximum gain)
		a = new float[order+1];
		a[0] = w1 / (float)Math.PI;
	    for (int n=1; n <= m; n++)
		    a[n] = (float)Math.sin(n*w1)*(float)Math.cos(n*w0)*win[n]/(n*(float)Math.PI);
		// shift impulse response to make filter causal:
		for (int n=m+1; n<=order; n++) a[n] = a[n - m];
		for (int n=0; n<=m-1; n++) a[n] = a[order - n];
		a[m] = w1 / (float)Math.PI;
	}

	public float[] filterGain () {

		// filter gain at uniform frequency intervals
		float[] g = new float[freqPoints+1];
		float theta, s, c, sac, sas;
		float gMax = -100.0f;
		float sc = 10.0f/(float)Math.log(10.0f);
		float t = (float)Math.PI / freqPoints;
		for (int i=0; i<=freqPoints; i++) {
			theta = i*t;
			sac = 0.0f;
			sas = 0.0f;
			for (int k=0; k<=order; k++) {
				c = (float)Math.cos(k*theta);
				s = (float)Math.sin(k*theta);
				sac += c*a[k];
				sas += s*a[k];
			}
			g[i] = sc*(float)Math.log(sqr(sac) + sqr(sas));
			gMax = Math.max(gMax, g[i]);
		}
		// normalise to 0 dB maximum gain
		for (int i=0; i<=freqPoints; i++) g[i] -= gMax;
        // normalise coefficients
        float normFactor = (float)Math.pow(10.0, -0.05*gMax);
        for (int i=0; i<=order; i++) a[i] *= normFactor;
		return g;
	}

}