setNumerator(coeffs.b);
setDenominator(coeffs.a);
With that done, we can now make the usual other complementary functions - hpf, bpf, bsf for high-pass, bandpass and bandstop. I'm not going to go into the theory of making these prototype functions and pole/zero placement in this post, but will come on to that later. Let's keep with code for the moment.
High-pass (hpf)
Following a similar procedure to that for lpf, the function is realised as follows:
// fcf is the cutoff frequency in Hz public void hpf(double fcf, int order) { double nfcf = fcf / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf > 0 && nfcf <= 0.5); // greater than 0 and less than Nyquist List<Complex> poles = getPoles(order); // get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf) / Math.PI; ComplexPlane splane = new ComplexPlane(); splane.zeros = new List<Complex>(); splane.poles = new List<Complex>(new Complex[order]); for (int i = 0; i < order; i++) { splane.poles[i] = w1 / poles[i]; splane.zeros.Add(new Complex(0, 0)); } ComplexPlane zplane = splane.BilinearTransform(); TransferFunction transferfunction = zplane.getTransferFunction(); Coefficients coeffs = transferfunction.getCoefficients(); double gain = transferfunction.getGainAt(new Complex(-1,0)); //gain at Nyquist // normalise the gain by dividing coeffs.b by gain coeffs.scaleGain(gain); setNumerator(coeffs.b); setDenominator(coeffs.a); }
This generates some good-looking sample graphs for different orders with the corner frequency again at 10kHz.
Bandpass (bpf)
This time we need to set two frequencies for the upper and lower frequencies
// fcf1 and fcf2 are the two cutoff frequencies in Hz. fcf1 < fcf2 public void bpf(double fcf1, double fcf2, int order) { double nfcf1 = fcf1 / PcmStreamSource.SampleRate; // normalised corner frequency double nfcf2 = fcf2 / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf1 > 0 && nfcf1 <= 0.5); // greater than 0 and less than Nyquist Debug.Assert(nfcf2 > 0 && nfcf2 <= 0.5); // greater than 0 and less than Nyquist Debug.Assert(fcf1 < fcf2); List<Complex> poles = getPoles(order); // get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf1) / Math.PI; double w2 = 2 * Math.PI * Math.Tan(Math.PI * nfcf2) / Math.PI; double w0 = Math.Sqrt(w1 * w2); double bw = w2 - w1; ComplexPlane splane = new ComplexPlane(); splane.poles = new List<Complex>(new Complex[order * 2]); splane.zeros = new List<Complex>(); for (int i = 0; i < order; i++) { Complex hba = poles[i] * (bw / 2.0); Complex temp = Complex.Sqrt(1 - Complex.Pow(w0 / hba, 2)); splane.poles[i] = hba * (1 + temp); splane.poles[order + i] = hba * (1 - temp); splane.zeros.Add(new Complex(0, 0)); } ComplexPlane zplane = splane.BilinearTransform(); TransferFunction transferfunction = zplane.getTransferFunction(); Coefficients coeffs = transferfunction.getCoefficients(); double centrefreq = (nfcf1 + nfcf2) / 2.0; Complex w = Complex.FromPolarCoordinates(1, 2 * Math.PI * centrefreq); double gain = transferfunction.getGainAt(w); // gain at centre freq // normalise the gain by dividing coeffs.b by gain coeffs.scaleGain(gain); setNumerator(coeffs.b); setDenominator(coeffs.a); }
Which looks like this for a bandpass with the lower frequency at 8kHz and the upper frequency at 12kHz.
Note, this has double the number of poles/zeros for the same order of filter compared to lpf or hpf.
Bandstop (bsf)
And finally a bandstop function:
// fcf1 and fcf2 are the two cutoff frequencies in Hz. fcf1 < fcf2 public void bsf(double fcf1, double fcf2, int order) { double nfcf1 = fcf1 / PcmStreamSource.SampleRate; // normalised corner frequency double nfcf2 = fcf2 / PcmStreamSource.SampleRate; // normalised corner frequency Debug.Assert(nfcf1 > 0 && nfcf1 <= 0.5); // greater than 0 and less than Nyquist Debug.Assert(nfcf2 > 0 && nfcf2 <= 0.5); // greater than 0 and less than Nyquist Debug.Assert(fcf1 < fcf2); List<Complex> poles = getPoles(order); // get warped freq double w1 = 2 * Math.PI * Math.Tan(Math.PI * nfcf1) / Math.PI; double w2 = 2 * Math.PI * Math.Tan(Math.PI * nfcf2) / Math.PI; double w0 = Math.Sqrt(w1 * w2); double bw = w2 - w1; ComplexPlane splane = new ComplexPlane(); splane.poles = new List<Complex>(new Complex[order * 2]); splane.zeros = new List<Complex>(new Complex[order * 2]); for (int i = 0; i < order; i++) { Complex hba = (bw / 2.0) / poles[i]; Complex temp = Complex.Sqrt(1 - Complex.Pow(w0 / hba, 2)); splane.poles[i] = hba * (1 + temp); splane.poles[order + i] = hba * (1 - temp); // 2n zeros at +-w0 splane.zeros[i] = new Complex(0, w0); splane.zeros[order + i] = new Complex(0, -w0); } ComplexPlane zplane = splane.BilinearTransform(); TransferFunction transferfunction = zplane.getTransferFunction(); Coefficients coeffs = transferfunction.getCoefficients(); double dcgain = transferfunction.getGainAt(Complex.One); // gain at dc double hfgain = transferfunction.getGainAt(new Complex(-1,0)); // gain at hf double gain = Math.Sqrt(dcgain * hfgain); // normalise the gain by dividing coeffs.b by gain coeffs.scaleGain(gain); setNumerator(coeffs.b); setDenominator(coeffs.a); }
Which looks like this:
Note, this has double the number of poles/zeros for the same order of filter compared to lpf or hpf.
That's a full complement of useful filter functions with the Butterworth characteristic.
No comments:
Post a Comment