java频谱补零的练习

1，生成一段模拟信号，进行模拟采样

/**
     * 生成一段模拟信号，模拟采样
     * @param fs 采样频率
     * @param n 采样点数
     * @param wantN 希望参与傅里叶变换的长度，因为想测试补零，所以该值会大于n
     * 模拟的这段信号包括的频率有10hz, 15hz, 20hz, 26.5hz，这段信号包括直流信号3
     * 注意：根据采样定理，采样频率必须大于信号频率的2倍
     * 因此，传递过来的采样频率务必要大于26.5 * 2
     */
    private double[] generateWaveFillZero(float fs, int n, int wantN) {
        double pi = Math.PI;
        double[] res = new double[wantN];
        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));
        }
        //后面补零
        for(int i = n; i < wantN; i++){
            res[i] = 0;
        }
        return res;
    }

2，对采集信号生成频谱

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

        //对采样结果进行傅里叶变换，使用commons-math包
        FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
        Complex[] result = fft.transform(wave, TransformType.FORWARD);
        log(Arrays.toString(result));

        //注意：不管是否补零，计算振幅的时候都是除以采样点数，而不是参与傅里叶变换的点数
        //计算第一个傅里叶变换结果的振幅，
        Complex complex0 = result[0];
        double amp0 = Math.sqrt(Math.abs(complex0.getReal()) * Math.abs(complex0.getReal())
                + Math.abs(complex0.getImaginary()) * Math.abs(complex0.getImaginary())) * 1.0 / N;//注意这里，哪怕是补零，这里还是除以N
        log(amp0);

        //计算剩余傅里叶变换结果的振幅,正式应用中只用计算到N / 2即可，因为是堆成
        for(int i = 1; i < WANT_N; i++){
            Complex complex = result[i];
            double amp = Math.sqrt(Math.abs(complex.getReal()) * Math.abs(complex.getReal())
                    + Math.abs(complex.getImaginary()) * Math.abs(complex.getImaginary())) * 2.0 / N;//注意这里，哪怕是补零，这里还是除以N
            log(amp);
        }
    }

 

3，用excel画图

分析：频率分辨率 = 100 / 1024；第272个点的频率 = （272 - 1）* 100 / 1024 = 26.46hz,正好对应的是4 * Math.sin(2 * pi * 26.5 * t + Math.toRadians(60)，振幅也正好能对应上。其他类推，比如第206个点的频率 = （206 - 1）* 100 /1024 = 20.02hz，正好对应的是3 * Math.cos(2 * pi * 20 * t + Math.toRadians(-30))，振幅也能对应上。