java频谱计算的工具类

直接上代码，逻辑还算比较清楚

public class SpectrumHelper {

    /**
     * 计算频谱
     * @param data 采样数据
     * @param samplingInterval 采样间隔,以秒为单位
     * @param needFillZero 是否需要补零
     * @param flyLen 参与傅里叶变换的长度
     */
    public static SpectrumBean caculateSpectrum(double[] data, double samplingInterval, boolean needFillZero, int flyLen){
        if(data == null || data.length == 0){
            return null;
        }
        int dataLen = data.length;
        double[] partakeFLYArr = null;//参与傅里叶变换的数组
        if(needFillZero){ //需要补零
            if(flyLen < dataLen){
                return null;
            }
            partakeFLYArr = new double[flyLen];
            System.arraycopy(data, 0, partakeFLYArr, 0, dataLen);
            for(int i = dataLen; i < flyLen; i++){
                partakeFLYArr[i] = 0;
            }
        }else{
            partakeFLYArr = data;
            flyLen = dataLen;
        }
        //进行傅里叶变换,使用的是commons-math包里面的快速傅里叶算法
        FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
        Complex[] complexArr = fft.transform(partakeFLYArr, TransformType.FORWARD);

        SpectrumBean spectrumBean = new SpectrumBean();
        double[] amplitude = new double[flyLen / 2];

        //计算振幅,除以的都是采样点的长度,而不是傅里叶变化的长度
        Complex complex = complexArr[0];
        amplitude[0] = Math.sqrt(Math.abs(complex.getReal()) * Math.abs(complex.getReal())
                + Math.abs(complex.getImaginary()) * Math.abs(complex.getImaginary())) * 1.0 / dataLen;
        for(int i = 1; i < flyLen / 2; i++){
            complex = complexArr[i];
            amplitude[i] = Math.sqrt(Math.abs(complex.getReal()) * Math.abs(complex.getReal())
                    + Math.abs(complex.getImaginary()) * Math.abs(complex.getImaginary())) * 2.0 / dataLen;
        }
        spectrumBean.setAmplitude(amplitude);//设置振幅

        //计算频率分辨率，除以的是傅里叶变换的长度，而不是采样点的长度
        double Fs = 1.0 / samplingInterval;//采样频率
        spectrumBean.setFrequencyResolution(Fs / flyLen);//设置频率分辨率

        return spectrumBean;
    }

}

SpectrumBean类表示的是频谱计算的结果，类定义如下：

/**
 * 用于保存频谱计算的结果
 */
public class SpectrumBean {

    /**
     * 频率分辨率
     */
    private double frequencyResolution;

    /**
     * 保存振幅信息,长度是参与傅里叶变换的长度,不是采样点的长度
     */
    private double[] amplitude;

    public double getFrequencyResolution() {
        return frequencyResolution;
    }

    public void setFrequencyResolution(double frequencyResolution) {
        this.frequencyResolution = frequencyResolution;
    }

    public double[] getAmplitude() {
        return amplitude;
    }

    public void setAmplitude(double[] amplitude) {
        this.amplitude = amplitude;
    }
}

测试的代码如下：

   /**
     * 正常的输出波形，正常的傅里叶变换，没有补零的操作
     * 频率分辨率 = 采样频率 / 参与傅里叶变换的点数，注意是参与傅里叶变换的点数而非采样点数，因为可能补零
     * 傅里叶变换结果的频率计算，第n个点的频率 = 频率分辨率 * （n - 1）
     * 所谓频谱，横轴是频率，纵轴是振幅,也就是频率振幅谱
     */
    @Test
    public void test3(){
        //模拟采样
        int Fs = 100;//采样频率
        int N = 256;//采样点数
        double[] wave = generateWaveNormal(Fs, N);
        log(Arrays.toString(wave));

        //对采样结果做频谱分析
        SpectrumBean spectrumBean = SpectrumHelper.caculateSpectrum(wave, 1.0f / Fs, false, 0);
        /*for(int i = 0; i < spectrumBean.getAmplitude().length; i++){
            log(spectrumBean.getAmplitude()[i]);
        }*/
        log(spectrumBean.getFrequencyResolution());

    }


    /**
     * 后面补零
     * 频率分辨率 = 采样频率 / 参与傅里叶变换的点数，注意是参与傅里叶变换的点数而非采样点数，
     * 在补零的情况下，参与傅里叶变换的点数是大于采样点数
     * 傅里叶变换结果的频率计算，第n个点的频率 = 频率分辨率 * （n - 1）
     * 所谓频谱，横轴是频率，纵轴是振幅,也就是频率振幅谱
     */
    @Test
    public void test4(){
        //模拟采样
        int Fs = 100;//采样频率
        int N = 256;//采样点数
        int WANT_N = 1024;//参与傅里叶变换的点数，实际上只采了256个点，后面的都补零了
        double[] wave = generateWaveNormal(Fs, N);
        log(Arrays.toString(wave));

        //对采样结果做频谱分析
        SpectrumBean spectrumBean = SpectrumHelper.caculateSpectrum(wave, 1.0 / Fs, true, WANT_N);
        /*for(int i = 0; i < spectrumBean.getAmplitude().length; i++){
            log(spectrumBean.getAmplitude()[i]);
        }*/
        log(spectrumBean.getFrequencyResolution());

    }

    
    private void log(Object obj){
        System.out.println(obj);
    }

    
    /**
     * 生成一段模拟信号，模拟采样
     * @param fs 采样频率
     * @param n 采样点数
     * 模拟的这段信号包括的频率有10hz, 15hz, 20hz, 26.5hz，这段信号包括直流信号
     * 注意：根据采样定理，采样频率必须大于信号频率的2倍
     * 因此，传递过来的采样频率务必要大于26.5 * 2
     */
    private double[] generateWaveNormal(float fs, int n) {
        double pi = Math.PI;
        double[] res = new double[n];
        float interval = 1 / fs;//采样间隔，单位是s
        for (int i = 0; i < n; i++) {
            double t = interval * i;
            res[i] = 3 + 1 * Math.cos(2 * pi * 10 * t) + 2 * Math.sin(2 * pi * 15 * t + Math.toRadians(30))
                    + 3 * Math.cos(2 * pi * 20 * t + Math.toRadians(-30)) + 4 * Math.sin(2 * pi * 26.5 * t + Math.toRadians(60));
        }
        return res;
    }

测试的结果打印到控制台，然后复制到excel种，使用excel画图

不补零的图如下：

补零的图如下：

结果分析

横坐标频率计算，比如第n个点的频率 = (n -1) * 频率分辨率，频率分辨率的计算见上面

纵坐标表示的是振幅。

频谱图上的信号和原始信号是能正好对应上的。