Java编程EMD

IMF需要满足的俩个条件：
1、极值点和过零点的数目应相等，或最多差一个
2、局部最大值和局部最小值的上下包络线均值为0

EMD的步骤：
1、包络线的获取：确定原始信号的极大值和极小值序列，采用三次样条曲线对极值点进行拟合，形成上下包络线；
注意用希尔伯特变换求取包络线的流程：

2、残余信号获取：原始信号x（t），平均信号e（t）做差，得到残余信号；

3、IMF条件检验，是否满足上面的2个条件

4、残差信号获取：原信号与平均信号作差，得到残余信号，将残差信号作为原始信号重复上面过程，不断得到各个固有模态函数：

5、分解中值：当残差信号为单调函数时，终止分解，此时原始信号可以用分解出来的iMF进行重构。

希尔伯特流程：

Hilbert 方法过程：

一个带求包络的信号x(t)，进行Hilbert变换的好 HHT(x(t)), 合成一个信号 x(t) + j*HHT( x(t) ), 然后对这个合成的信号取幅值部分 y(t)= Amp( [x(t) + j*HHT( x(t) )] ), 此时y(t)就是得到的上包络，下包络y'(t) = -y(t).

 针对上面的流程进行编程：

package model.EXETest;
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class HBTTools {

        // 计算N点DFT的一半，用于递归实现完整DFT
        public static void dftRecursive(double[] real, double[] imag, int N, double[] omega) {
            if (N == 1) return;

            // 递归调用计算一半的DFT
            dftRecursive(real, imag, N / 2, omega);

            // 对称复制剩余的一半并递归调用
            for (int k = 0; k < N / 2; k++) {
                real[k + N / 2] = real[k];
                imag[k + N / 2] = -imag[k];
            }

            // 递归过程中进行按频率的循环位移
            double tmpReal, tmpImag;
            for (int k = 0; k < N / 2; k++) {
                double kth = omega[N / 2] * k;
                tmpReal = real[N / 2 + k];
                tmpImag = imag[N / 2 + k];
                real[N / 2 + k] = real[k] + tmpReal * Math.cos(kth) - tmpImag * Math.sin(kth);
                imag[N / 2 + k] = imag[k] + tmpReal * Math.sin(kth) + tmpImag * Math.cos(kth);

                real[k] = real[k] + tmpReal * Math.cos(kth) + tmpImag * Math.sin(kth);
                imag[k] = imag[k] - tmpReal * Math.sin(kth) + tmpImag * Math.cos(kth);
            }
        }

        public static double[] hartleyTransform(double[] x) {
            int N = x.length;
            double[] real = new double[N];
            double[] imag = new double[N];
            double[] omega = new double[N];

            // 初始化omega数组
            for (int n = 0; n < N; n++) {
                omega[n] = -2 * Math.PI * n / N;
            }

            // 复数部分初始化为0
            Arrays.fill(imag, 0.0);

            // 实部赋值为输入序列
            System.arraycopy(x, 0, real, 0, N);

            // 递归计算DFT
            dftRecursive(real, imag, N, omega);

            // 输出DFT结果
            for (int n = 0; n < N; n++) {
                x[n] = Math.sqrt(real[n] * real[n] + imag[n] * imag[n]); // 计算幅度
            }
            return x;
        }


    //去取文件获取数据返回集合Double
    private static List<Double> readFile(String path) {
        List<Double> ll=new ArrayList<>();
        String line="";
        try
        {
            BufferedReader in=new BufferedReader(new FileReader(path));
            while ((line=in.readLine()) != null)
            {
                //表示字符串长度不受限制
//                String[] str=line.split(" ",-1);
                //先把外面的替换掉
                String c = line.replaceAll("\r|\n", "");
                String[] str=c.split(" ");
                for(int i=0;i<str.length;i++){
                    if (i % 5 == 0 || i % 5 == 1 || i % 5 == 2) {
                        double dat = Double.parseDouble(str[i]);
                        ll.add(dat);
                    }
                }
            }
            in.close();
        } catch (IOException e) {
            e.printStackTrace();
        }

        return ll;
    }

    //信号取反，获取下包络
    public static double[] getDownLineData(double[] dek){
        double[] buk=new double[dek.length];
        for(int i=0;i<dek.length;i++){
            buk[i] = dek[i]*(-1);
        }
        return dek;
    }

    //信号取平均 h(t)=x(t)-e(t);
    public static double[] diffentenceData(double[] orangin,double[] up,double[] dowm){
        double[] result=new double[up.length];
        for(int j=0;j<up.length;j++){
            double value=(up[j]-dowm[j])/2;
            result[j]=orangin[j]-value;
        }
        return result;
    }
    //结果获取
    public static double[] IMFR(double[] orginle){
        // 希尔伯特变换 获取上包络
        double[] up=hartleyTransform(orginle);

        //对数据取反，获取下包络
        double[] backOR=getDownLineData(orginle);
        //获取下包络
        double[] dowmN=hartleyTransform(backOR);
        //再取反，获取原信号下包络
        double[] down=getDownLineData(dowmN);

        //原信号减去包络的均值，获取IMF信号
        double[] imf=diffentenceData(orginle,up,down);
        return imf;
    }
        public static void main(String[] args) {
            // 示例序列
            double[] x = {1, 2, 3, 4, 5, 6, 7, 8,9,8,7,6,5,4,3,2,1};
            String path="C:\\Users\\seer\\Desktop\\测量无效\\sport\\6327894be0924b44ab82a45c737ee932_1713413035_1713413225963_160_1.txt";
            List<Double> fileData=readFile(path);
            double[] base=SystemTools.getList(fileData);

            //传入原始信号
            double[] imf=IMFR(base);
            for(int k=0;k<8;k++){
                imf=IMFR(imf);
                // 输出变换后的序列
            }
            StringBuilder sb=new StringBuilder();
            for(int i=0;i<imf.length-1;i++){
                sb.append(imf[i]+" ");
            }
            System.out.print(sb.toString());

        }

}