using CommonLib.LOG; using NWaves.FeatureExtractors; using NWaves.FeatureExtractors.Multi; using NWaves.FeatureExtractors.Options; using System; using System.Data.SqlTypes; using System.IO; using System.Linq; using System.Text; using System.Windows.Forms; namespace DataCollectionSystem.CommonLib.Data { public class DataProcess { alglib.clusterizerstate s1;//释放聚类空间 alglib.kmeansreport rep1;//聚类结果报告 alglib.clusterizerstate s2; alglib.kmeansreport rep2; public int repeat = 0;//次数 //定义固定窗 int MatrixLen = 20; //数据窗的长度(行) int MaxtrixDimension =4;//数据的特征维度(列) double[,] fixed_window = null;//new double[MatrixLen, MaxtrixDimension]; //定义特征展示窗 double[] features = null;//new double[35]; //int Len = 3600; //回溯1小时Mel int Len = 300; //回溯1小时Mel public double[] Mel1 = null;//new double[Len];//特征窗 public double[] Mel2 = null;//new double[Len];//特征窗 public double[] Time1 = null;//new double[Len];//特征窗 public double[] Time2 = null;//new double[Len];//特征窗 public double[] spectral1 = null;//new double[Len];//特征窗 public double[] spectral2 = null;//new double[Len];//特征窗 public double[,] all_features = new double[35, 300];//new double[Len];//特征窗 public int pos = 0; //定义一个数据窗 int DataLen = 20; //数据窗的长度(行) int DataDimension = 4;//数据的特征维度(列) //double[,] Data_window = null;// new double[DataLen, DataDimension]; double[,] Data_window = null;//new double[DataLen, DataDimension]; //固定莱特窗 int WrightWindowLen =20; double[] Wright = null;//new double[WrightWindowLen];//莱特窗 double k = 8;//倍数 //定义窗的指标 double avr = 0; //int sum = 0; //double sum1 = 0; double Standard = 0; //double ss = 0; double Residual = 0; /// /// 把环形缓冲区的数据计算出来 /// /// /// public Boolean checkError(string identName,double[] RingArray1, ref Boolean anomaly, ref int anomalyInt )//double Feature_0, double Feature_1) { int len = RingArray1.Length; float[] RingArray = new float[len]; for (int i = 0; i < len; i++) { RingArray[i] = (float)RingArray1[i]; } Boolean result = true; alglib.clusterizerstate s1;//释放聚类空间 alglib.kmeansreport rep1;//聚类结果报告 alglib.clusterizerstate s2; alglib.kmeansreport rep2; if (Data_window == null) { //定义固定窗 fixed_window = new double[MatrixLen, MaxtrixDimension]; //定义特征展示窗 features = new double[35]; Mel1 = new double[Len];//特征窗 Mel2 = new double[Len];//特征窗 Time1 = new double[Len];//特征窗 Time2 = new double[Len];//特征窗 spectral1 = new double[Len];//特征窗 spectral2 = new double[Len];//特征窗 //定义一个数据窗 Data_window = new double[DataLen, DataDimension]; //固定莱特窗 Wright = new double[WrightWindowLen];//莱特窗 } #region 特征提取 var mfccOptions = new MfccOptions //梅尔倒谱系数提取 { SamplingRate = 32000, FeatureCount = 13, FrameDuration = 1, HopDuration = 1, FilterBankSize = 26, // PreEmphasis = 0.97, }; var mfccExtractor = new MfccExtractor(mfccOptions); //var mfccVectors = mfccExtractor.ComputeFrom(fArr); //得到mel系数 var mfccVectors = (mfccExtractor.ComputeFrom(RingArray))[0]; //得到mel系数 //float[] mfccvectors = mfccVectors[0]; //保存mel特征 for (int i = 0; i < 13; ++i) { features[i] = mfccVectors[i]; //fArr[i] = f; } var lpccOptions = new LpccOptions //LPCC系数提取 { SamplingRate = 32000, FeatureCount = 10, FrameDuration = 1, HopDuration = 1, LpcOrder = 1, //LifterSize = _lifterSize, // PreEmphasis = _preEmphasis, //Window = _window }; var LpccExtractor = new LpccExtractor(lpccOptions); var LpccVectors = (LpccExtractor.ComputeFrom(RingArray))[0]; //得到前10维LPCC系数 //float[] Lpccvectors = LpccVectors[0]; for (int i = 0; i < 10; ++i) { features[i + 13] = LpccVectors[i]; //fArr[i] = f; } var spectralOptions = new MultiFeatureOptions //频域特征提取 { SamplingRate = 32000, //FeatureList = spectralFeatureSet, FrameDuration = 1, HopDuration = 1, //FftSize = _blockSize, //Frequencies = _frequencies, //PreEmphasis = _preEmphasis, //Window = _window, //Parameters = _parameters }; var Spectral_featureExtractor = new SpectralFeaturesExtractor(spectralOptions); var Spectral_featureVectors = (Spectral_featureExtractor.ComputeFrom(RingArray))[0]; //频域特征向量 for (int i = 0; i < 8; i++) { features[i + 23] = Spectral_featureVectors[i]; //fArr[i] = f; } var timeOptions = new MultiFeatureOptions //时域特征提取 { SamplingRate = 31998, FrameDuration = 1, HopDuration = 1, //FeatureList = string.Join(",", FeatureDescriptions), //Parameters = _parameters }; var Time_featureExtractor = new TimeDomainFeaturesExtractor(timeOptions); //float[] Time_HardArr = new float[10001]; //for (int i = 0; i < 10000; i++) //{ //Time_HardArr[i] = HardArr[i]; //fArr[i] = f; //} //Time_HardArr[10000] = HardArr[9999]; var timeVectors = (Time_featureExtractor.ComputeFrom(RingArray))[0]; //时域特征向量 //var timeVectors = Time_featureExtractor.ComputeFrom(HardArr); //时域特征向量 for (int i = 0; i < 4; ++i) { features[i + 31] = timeVectors[i]; //fArr[i] = f; } //all //int Len = 300; //回溯1小时Mel //features = new double[35]; // public double[,] all_features = new double[35, 300];//new double[Len];//特征窗 #endregion #region 特征窗更新 double m1 = features[0]; for (int i = 1; i < Len; i++)//窗更新 { Mel1[i - 1] = Mel1[i]; } Mel1[Len - 1] = m1; double m2 = features[15]; for (int i = 1; i < Len; i++)//窗更新 { Mel2[i - 1] = Mel2[i]; } Mel2[Len - 1] = m2; double t1 = features[31]; for (int i = 1; i < Len; i++)//窗更新 { Time1[i - 1] = Time1[i]; } Time1[Len - 1] = t1; double t2 = features[21]; for (int i = 1; i < Len; i++)//窗更新 { Time2[i - 1] = Time2[i]; } Time2[Len - 1] = t2; double sp1 = features[23]; for (int i = 1; i < Len; i++)//窗更新 { spectral1[i - 1] = spectral1[i]; } spectral1[Len - 1] = sp1; double sp2 = features[25]; for (int i = 1; i < Len; i++)//窗更新 { spectral2[i - 1] = spectral2[i]; } spectral2[Len - 1] = sp2; #endregion // 选取特征 阈值a b c d e double a = features[31]; double b = features[32]; double c = features[33]; //double d = features[23];//对质心不用归一 double d = features[23] / 32000;//对质心归一 double e = features[29]; try { //方法1 阈值 //if(repeat > 1 ) //{ //string Info1 = $"查看特征:由于短时能量:{a},RMS:{b},过零率:{c},质心:{d},谱熵:{e}"; //LogStorage.Logs.Add(LogLevel.Error, Info1); //} // //if (a > 0.015 || b > 0.12 || c > 0.14 || d > 0.3) if (a > 1.5 && b > 1.2 && c > 1.4 && d > 3) //if (a > 0.015 && b > 0.12 && c > 0.14 && d > 1600 )//残差减去5倍标准差 //if (a > 0.15 || b > 0.12 || c > 0.14 || d > 1600 )//残差减去5倍标准差 //if(j>3&&laiTe[j-1]> laiTe[j]&& laiTe[j - 1]> laiTe[j - 2]&& ((laiTe_2[j-1]/ laiTe[j - 1])>0.1)) { /* if(a > 0.015) { string Info1 = $"法1(阈值判别):由于C大于设定的阈值(0.015),在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error, Info1); } if (b > 0.12) { string Info1 = $"法1(阈值判别):由于RMS:{b}大于设定的阈值(0.12),在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error, Info1); } if (c > 0.14) { string Info1 = $"法1(阈值判别):由于过零率:{c}大于设定的阈值(0.14),在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error, Info1); } if (d > 0.3) { string Info1 = $"法1(阈值判别):由于谱熵:{e}大于设定的阈值(0.3),在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error, Info1); }*/ /* Console.ForegroundColor = ConsoleColor.Red; Console.WriteLine("******************************************************异常出现******************************************************"); Console.WriteLine("法1(阈值判别):由于短时能量:{0},RMS:{1},过零率:{2},质心:{3},谱熵:{4} 之中有大于设定的阈值,在第{5}点出现了疑似异常", a, b, c, d, e, j); Console.WriteLine("******************************************************异常出现******************************************************"); */ string Info1 = $"{identName} 法1(阈值判别):由于短时能量:{a}大于设定的阈值(0.015),在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error1, Info1); anomaly = true; anomalyInt = 1; //string Info1 = $"法1(阈值判别):由于短时能量:{a},RMS:{b},过零率:{c},质心:{d},谱熵:{e} 之中有大于设定的阈值,在第{repeat}点出现了疑似异常"; //LogStorage.Logs.Add(LogLevel.Error, Info1); } //Console.WriteLine("短时能量:{0},RMS:{1},过零率:{2},质心:{3},谱熵:{4} ,目前运行到了第{5}点 运行正常!!!", a, b, c, d, e,j); //Console.WriteLine("目前运行到了第{0}点", j); //Console.ForegroundColor = ConsoleColor.White; //Console.ReadKey(); //方法2 均值聚类中心差异判别 //固定窗 fixed_window //int MatrixLen = 20; //数据窗的长度(行) //int MaxtrixDimension = 4;//数据的特征维度(列) //double[,] fixed_window = new double[MatrixLen, MaxtrixDimension]; for (int z = 1; z < MatrixLen; z++)//窗更新 { fixed_window[z - 1, 0] = fixed_window[z, 0]; fixed_window[z - 1, 1] = fixed_window[z, 1]; fixed_window[z - 1, 2] = fixed_window[z, 2]; fixed_window[z - 1, 3] = fixed_window[z, 3]; } //Wright[Wright.Length - 1] = E; fixed_window[MatrixLen - 1, 0] = a; fixed_window[MatrixLen - 1, 1] = b; fixed_window[MatrixLen - 1, 2] = c; fixed_window[MatrixLen - 1, 3] = d; if (repeat >= 20 && repeat % 20 == 0) //满足取窗条件 { //开始聚类 alglib.clusterizercreate(out s2); alglib.clusterizersetpoints(s2, fixed_window, 2); alglib.clusterizersetkmeanslimits(s2, 2000, 0); alglib.clusterizerrunkmeans(s2, 2, out rep2);//聚2类 double a1 = rep2.c[0, 0]; double a2 = rep2.c[0, 1]; double a3 = rep2.c[0, 2]; double a4 = rep2.c[0, 3]; double b1 = rep2.c[1, 0]; double b2 = rep2.c[1, 1]; double b3 = rep2.c[1, 2]; double b4 = rep2.c[1, 3]; double U1 = Math.Abs(a1 - b1); // + Math.Abs(a2 - b2) + Math.Abs(a3 - b3); double U2 = Math.Abs(a2 - b2); double U3 = Math.Abs(a3 - b3); double U4 = Math.Abs(a4 - b4); int sum = 0; for (int i = 0; i < rep2.cidx.Length; i++) { //Console.WriteLine(Wright[i]); sum += rep2.cidx[i]; } // int sum = rep2.cidx.Sum(); if (U1 > 1 && U2 > 5 && U3 > 7 && U4 > 1.2 && sum > 9) //if ((U1 > 0.01 || U2 > 0.05 || U3 > 0.07 || U4 > 0.12) && sum > 9) //if (U1 > 0.01 && U2 > 0.05 && U3 > 0.07 && U4 > 800 && sum > 9)//残差减去5倍标准差 //if (U1 > 0.01 || U2 > 0.05 || U3 > 0.07 || U4 > 800)//残差减去5倍标准差 //if (a > 0.15 || b > 0.12 || c > 0.14 || d > 1600 )//残差减去5倍标准差 //if(j>3&&laiTe[j-1]> laiTe[j]&& laiTe[j - 1]> laiTe[j - 2]&& ((laiTe_2[j-1]/ laiTe[j - 1])>0.1)) { /* Console.ForegroundColor = ConsoleColor.Blue; Console.WriteLine("******************************************************异常出现******************************************************"); Console.WriteLine("法2(均值聚类):由于聚类中心c1({0},{1},{2},{3})与聚类中心c2({4},{5},{6},{7})相差过大,在第{8}点出现了疑似异常", a1, a2, a3, a4, b1, b2, b3, b4, j); Console.WriteLine("第{0}聚类,异常的个数为{1}", j, sum); Console.WriteLine("******************************************************异常出现******************************************************"); */ string Info1 = $"{identName} 法2(均值聚类):由于聚类中心c1({a1},{a2},{a3},{a4})与聚类中心c2({b1},{b2},{b3},{b4})相差过大,在第{repeat}点出现了疑似异常,异常个数{sum}"; LogStorage.Logs.Add(LogLevel.Error2, Info1); anomaly = true; anomalyInt = 1; } //Console.WriteLine("短时能量:{0},RMS:{1},过零率:{2},质心:{3},谱熵:{4} ,目前运行到了第{5}点 运行正常!!!", a, b, c, d, e, j); //Console.ForegroundColor = ConsoleColor.White; //Console.ReadKey(); } // 滑动窗口判别 //对数据窗更新 //int DataLen = 20; //数据窗的长度(行) //int DataDimension = 4;//数据的特征维度(列) //double[,] Data_window = new double[DataLen, DataDimension]; for (int z = 1; z < DataLen; z++)//窗更新 { Data_window[z - 1, 0] = Data_window[z, 0]; Data_window[z - 1, 1] = Data_window[z, 1]; Data_window[z - 1, 2] = Data_window[z, 2]; Data_window[z - 1, 3] = Data_window[z, 3]; } //Wright[Wright.Length - 1] = E; Data_window[DataLen - 1, 0] = a; Data_window[DataLen - 1, 1] = b; Data_window[DataLen - 1, 2] = c; Data_window[DataLen - 1, 3] = d; alglib.clusterizercreate(out s1); alglib.clusterizersetpoints(s1, Data_window, 2); alglib.clusterizersetkmeanslimits(s1, 10, 0); alglib.clusterizerrunkmeans(s1, 1, out rep1);//聚1类 alglib.clusterizercreate(out s2); alglib.clusterizersetpoints(s2, Data_window, 2); alglib.clusterizersetkmeanslimits(s2, 2000, 0); alglib.clusterizerrunkmeans(s2, 2, out rep2);//聚2类 double E = Math.Abs(rep1.energy - rep2.energy); //double E = Math.Abs(rep1.energy - rep2.energy);//只检测突然增大 for (int i = 1; i < Wright.Length; i++)//窗更新 { Wright[i - 1] = Wright[i]; } Wright[Wright.Length - 1] = E; double sum1 = 0; for (int i = 0; i < Wright.Length; i++) { //Console.WriteLine(Wright[i]); sum1 += Wright[i]; } //sum1 = Wright.Sum(); avr = sum1 / Wright.Length;//均值 double ss = 0; for (int i = 0; i < Wright.Length; i++) { ss += Math.Pow((Wright[i] - avr), 2); } //biaoZhun = Math.Sqrt(s/(laiTe.Length-1));//除以N-1 Standard = Math.Sqrt(ss / (Wright.Length));//除以N Residual = Math.Abs(Wright[Wright.Length - 1] - avr); if (Residual > k * Standard && a > 1 && repeat > 40) {/* Console.ForegroundColor = ConsoleColor.Yellow; Console.WriteLine("******************************************************异常出现******************************************************"); Console.WriteLine("法3(自适应聚类):由于残差{0}大于标准差{1},在第{2}点出现了疑似异常", Residual, Standard, j); //Console.WriteLine("第{0}聚类,异常的个数为{1}", j, sum); Console.WriteLine("******************************************************异常出现******************************************************"); */ string Info1 = $"{identName} 法3(自适应聚类):由于残差{Residual}大于{k}倍标准差{Standard},在第{repeat}点出现了疑似异常"; LogStorage.Logs.Add(LogLevel.Error, Info1); anomaly = true; anomalyInt = 1; /* Info1 = $"第{j}聚类,异常的个数为{sum}"; LogStorage.Logs.Add(LogLevel.Error, Info1); */ } //if(repeat>20) //{ //string Info1 = $"查看计算值:和1{sum1}、和2{ss}、平均值{avr}、第{repeat}点"; //LogStorage.Logs.Add(LogLevel.Error, Info1); //} repeat++; if (repeat > 8640000) { repeat = 0; } //保存特征值 if(pos<300) { // all_features features int vector_size = features.Count(); for (int i = 0; i < vector_size; i++) { all_features[i, pos] = features[i]; } pos++; } else { if (pos == 300) { SaveCSV(all_features); pos++; } } } catch (Exception ex) { LogStorage.Logs.Add(LogLevel.Error, ex.Message + ex.StackTrace); } // } return result; } private void SaveCSV(double[,] all_features) { int len0 = all_features.GetLength(0); int len1 = all_features.GetLength(1); StringBuilder sb = new StringBuilder(); for(int i=0; i