Java代码实现 6. order butterworth filter with F(0) 20Hz.

下面是一个简单的 Java 代码实现，可以生成一个6阶的Butterworth滤波器，截止频率为20Hz。

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.function.Sigmoid;
import org.apache.commons.math3.analysis.function.Sinc;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.complex.ComplexUtils;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;
import java.util.ArrayList;

public class ButterworthFilter {
    private final int order;
    private final double cutoff;

    public ButterworthFilter(int order, double cutoff) {
        this.order = order;
        this.cutoff = cutoff;
    }

    public double[] apply(double[] input) {
        // 计算采样频率
        double sampleRate = 1.0 / input.length;

        // 计算截止频率在频域中的位置
        double normalizedCutoff = cutoff / (sampleRate / 2);

        // 构建巴特沃斯滤波器的极点
        ArrayList<Complex> poles = new ArrayList<>();
        for (int i = 0; i < order; i++) {
            double real = -Math.sin(Math.PI * (2 * i + 1) / (2 * order));
            double imag = Math.cos(Math.PI * (2 * i + 1) / (2 * order));
            poles.add(new Complex(real, imag));
        }

        // 构建巴特沃斯滤波器的传递函数
        UnivariateFunction transferFunction = new UnivariateFunction() {
            @Override
            public double value(double frequency) {
                Complex numerator = Complex.ONE;
                Complex denominator = Complex.ONE;

                for (Complex pole : poles) {
                    double real = pole.getReal();
                    double imag = pole.getImaginary();
                    Complex term = ComplexUtils.polarToCartesian(1, -2 * Math.PI * frequency * real, -2 * Math.PI * frequency * imag);
                    numerator = numerator.multiply(term);
                    denominator = denominator.multiply(new Complex(1, -2 * Math.PI * frequency * real, -2 * Math.PI * frequency * imag));
                }

                return numerator.divide(denominator).getReal();
            }
        };

        // 计算滤波器的频率响应
        int fftSize = input.length * 2;
        double[] frequencyResponse = new double[fftSize];
        for (int i = 0; i < fftSize; i++) {
            double frequency = i / (double)fftSize * sampleRate;
            frequencyResponse[i] = transferFunction.value(frequency);
        }

        // 将频率响应进行反变换，得到滤波器的时域响应
        FastFourierTransformer transformer = new FastFourierTransformer();
        double[] impulseResponse = transformer.transform(frequencyResponse, TransformType.INVERSE);

        // 将时域响应进行归一化，并截取前一半，得到滤波器的实际响应
        double[] filter = new double[input.length];
        for (int i = 0; i < input.length; i++) {
            filter[i] = impulseResponse[i] / input.length;
        }
        filter = java.util.Arrays.copyOfRange(filter, 0, input.length / 2);

        // 将滤波器应用到输入信号上
        double[] output = new double[input.length];
        for (int i = 0; i < input.length; i++) {
            double t = i * sampleRate;
            double sum = 0;
            for (int j = 0; j < filter.length; j++) {
                double tau = (j - filter.length / 2) * sampleRate;
                sum += filter[j] * new Sinc().value((t - tau) * normalizedCutoff);
            }
            output[i] = sum;
        }

        return output;
    }
}

在这个实现中，使用了 Apache Commons Math 库来进行一些数学计算，比如傅里叶变换。这个实现中的巴特沃斯滤波器是基于极点的实现方法，通过计算滤波器的传递函数，然后进行反傅里叶变换，得到滤波器的时域响应，最后将滤波器应用到输入信号上，得到滤波后的输出信号。