基于emd算法的实现算法.rar

/* 1.采样频率改为50Hz,即每秒钟采样50个点。 2.每40秒做一次EMD。 3.asp,hea为comp类型的结构体，分别存储着计算出来的呼吸和心跳的数据； 4.calculate()为计算呼吸率和心率的函数； 5.max_val函数，由于会检测出频率为0时的最大值，故从1开始往后比较，不确定是否准确； */ #include <stdio.h> #include <stdlib.h> #include <math.h> #define N 2048 #define PI 3.141592653 //采样频率为25.6Hz；如果包含两个周期的话，呼吸最长的周期大约八到九秒，所以512个点可以了； #define freq 50 #define peri 40.96 //以下是spline函数的一些变量的结构体声明； typedef enum _condition { NATURAL, CLAMPED, NOTaKNOT }condition; typedef struct _coefficient { double a3; double b2; double c1; double d0; }coefficient; typedef struct _point { coefficient *coe; double *xCoordinate; double *yCoordinate; double f0; double fn; int num; condition con; }point; typedef struct _comp // 复数结构体声明； { double re; double im; }comp; comp asp[N],hea[N]; // 存储的呼吸和心跳重构信号； /***************************极小值***************************/ int extr_min(double* y, int* b) { int j = 0; int i; //找到所有的极小值，并且将对应的横坐标存入数组中； for(i=1; i<N-1; i++) { if(y[i]<y[i+1]&&y[i]<y[i-1]) { b[++j] = i;//从第一个开始，第零个准备存原始信号的第一个数； } } return j; } /***************************极大值***************************/ int extr_max(double* y, int* b) { int j = 0; int i; //找到所有的极大值，并且将对应的横坐标存入数组中； for(i=1; i<N-1; i++) { if(y[i]>y[i+1]&&y[i]>y[i-1]) { b[++j] = i;//从第一个开始，第零个准备存原始信号的第一个数； } } return j; } /***************************FFT***************************/ void bitrev(comp data[N]) //fft的蝶形运算之前先进行排序； { int p=1, q, i; int bit_rev[ N ]; double xx_r[ N ]; bit_rev[ 0 ] = 0; while( p < N ) { for(q=0; q<p; q++) { bit_rev[ q ] = bit_rev[ q ] * 2; bit_rev[ q + p ] = bit_rev[ q ] + 1; } p *= 2; } for(i=0; i<N; i++) xx_r[ i ] = data[i].re; for(i=0; i<N; i++) data[i].re = xx_r[ bit_rev[i] ]; //xx_r[N]要定义全局变量； } void fft(comp data[N]) { int L; //碟形运算的每级 L = 1, 2,......,M int J; //J = 0, 1, 2....., 2^(L - 1) - 1 int B; //B = 2^(L - 1) int P; //P = J * 2^(M - L) int M = 11; // M的值取决于N，M = log2(N); int k; double datak_re,datak_im,datakb_re,datakb_im; bitrev(data); //排序； //碟形运算; for(L = 1; L <= M; L++) { B = 1 << (L - 1); for(J = 0; J <= B - 1; J++) { P = J * (1 << (M - L)); for(k = J; k <= N - 1; k += (1 << L)) { datak_re = data[k].re; datak_im = data[k].im; datakb_re = data[k+B].re; datakb_im = data[k+B].im; data[k].re = datak_re + cos(2*PI*P/N)*datakb_re + sin(2*PI*P/N)*datakb_im; data[k].im = datak_im + cos(2*PI*P/N)*datakb_im - sin(2*PI*P/N)*datakb_re; data[k+B].re = datak_re - cos(2*PI*P/N)*datakb_re - sin(2*PI*P/N)*datakb_im; data[k+B].im = datak_im - cos(2*PI*P/N)*datakb_im + sin(2*PI*P/N)*datakb_re; } } } } /**********************************spline********************************************/ int spline(point *point) { double *x = (*point).xCoordinate, *y = (*point).yCoordinate; int n = (*point).num - 1; coefficient *coe = (*point).coe; condition con = (*point).con; double *h, *d; double *a, *b, *c, *f, *m; double temp; int i; h = (double *)malloc( n * sizeof(double) ); /* 0,1--(n-1),n */ d = (double *)malloc( n * sizeof(double) ); /* 0,1--(n-1),n */ a = (double *)malloc( (n + 1) * sizeof(double) ); /* 特别使用,1--(n-1), n */ b = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1), n */ c = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1),特别使用 */ f = (double *)malloc( (n + 1) * sizeof(double) ); /* 0,1--(n-1), n */ m = b; if(f == NULL) { free(h); free(d); free(a); free(b); free(c); free(f); return 1; } /* 计算 h[] d[] */ for (i = 0; i < n; i++) { h[i] = x[i + 1] - x[i]; d[i] = (y[i + 1] - y[i]) / h[i]; /* printf("%f\t%f\n", h[i], d[i]); */ } /********************** ** 初始化系数增广矩阵 **********************/ a[0] = (*point).f0; c[n] = (*point).fn; /* 计算 a[] b[] d[] f[] */ switch(con) { case NATURAL: f[0] = a[0]; f[n] = c[n]; a[0] = 0; c[n] = 0; c[0] = 0; a[n] = 0; b[0] = 1; b[n] = 1; break; case CLAMPED: f[0] = 6 * (d[0] - a[0]); f[n] = 6 * (c[n] - d[n - 1]); a[0] = 0; c[n] = 0; c[0] = h[0]; a[n] = h[n - 1]; b[0] = 2 * h[0]; b[n] = 2 * h[n - 1]; break; case NOTaKNOT: f[0] = 0; f[n] = 0; a[0] = h[0]; c[n] = h[n - 1]; c[0] = -(h[0] + h[1]); a[n] = -(h[n - 2] + h[n - 1]); b[0] = h[1]; b[n] = h[n - 2]; break; } for (i = 1; i < n; i++) { a[i] = h[i - 1]; b[i] = 2 * (h[i - 1] + h[i]); c[i] = h[i]; f[i] = 6 * (d[i] - d[i - 1]); } /* for (i = 0; i < n+1; i++) printf("%f\n", c[i]); */ /************* ** 高斯消元 *************/ /* 第0行到第(n-3)行的消元 */ for(i = 0; i <= n - 3; i++) { /* 选主元 */ if( fabs(a[i + 1]) > fabs(b[i]) ) { temp = a[i + 1]; a[i + 1] = b[i]; b[i] = temp; temp = b[i + 1]; b[i + 1] = c[i]; c[i] = temp; temp = c[i + 1]; c[i + 1] = a[i]; a[i] = temp; temp = f[i + 1]; f[i + 1] = f[i]; f[i] = temp; } temp = a[i + 1] / b[i]; a[i + 1] = 0; b[i + 1] = b[i + 1] - temp * c[i]; c[i + 1] = c[i + 1] - temp * a[i]; f[i + 1] = f[i + 1] - temp * f[i]; } /* 最后3行的消元 */ /* 选主元 */ if( fabs(a[n - 1]) > fabs(b[n - 2]) ) { temp = a[n - 1]; a[n - 1] = b[n - 2]; b[n - 2] = temp; temp = b[n - 1]; b[n - 1] = c[n - 2]; c[n - 2] = temp; temp = c[n - 1]; c[n - 1] = a[n - 2]; a[n - 2] = temp; temp = f[n - 1]; f[n - 1] = f[n - 2]; f[n - 2] = temp; } /* 选主元 */ if( fabs(c[n]) > fabs(b[n - 2]) ) { temp = c[n]; c[n] = b[n - 2]; b[n - 2] = temp; temp = a[n]; a[n] = c[n - 2]; c[n - 2] = temp; temp = b[n]; b[n] = a[n - 2]; a[n - 2] = temp; temp = f[n]; f[n] = f[n - 2]; f[n - 2] = temp; } /* 第(n-1)行消元 */ temp = a[n - 1] / b[n - 2]; a[n - 1] = 0; b[n - 1] = b[n - 1] - temp * c[n - 2]; c[n - 1] = c[n - 1] - temp * a[n - 2]; f[n - 1] = f[n - 1] - temp * f[n - 2]; /* 第n行消元 */ temp = c[n] / b[n - 2]; c[n] = 0; a[n] = a[n] - temp * c[n - 2]; b[n] = b[n] - temp * a[n - 2]; f[n] = f[n] - temp * f[n - 2]; /* 选主元 */ if( fabs(a[n]) > fabs(b[n - 1]) ) { temp = a[n]; a[n] = b[n - 1]; b[n - 1] = temp; temp = b[n]; b[n] = c[n - 1]; c[n - 1] = temp; temp = f[n]; f[n] = f[n - 1]; f[n - 1] = temp; } /* 最后一次消元 */ temp = a[n] / b[n-1]; a[n] = 0; b[n] = b[n] - temp * c[n - 1]; f[n] = f[n] - temp * f[n - 1]; /* 回代求解 m[] */ m[n] = f[n] / b[n]; m[n - 1] = (f[n - 1] - c[n - 1] * m[n]) / b[n-1]; for(i = n - 2; i >= 0; i--) { m[i] = ( f[i] - (m[i + 2] * a[i] + m[i + 1] * c[i]) ) / b[i]; } /* for (i = 0; i < n+1; i++) printf("%f\n", m[i]); */ free(a); free(c); free(f); /************ ** 放置参数 ************/ for(i = 0; i < n; i++) { coe[i].a3 = (m[i + 1] - m[i]) / (6 * h[i]); coe[i].b2 = m[i] / 2; coe[i].c1 = d[i] - (h[i] / 6) * (2 * m[i] + m[i + 1]); coe[i].d0 = y[i]; } free(h); free(d); free(b); return n +