数组数据分析算法中峰区域的确定
做数据分析算法,使用MATLAB进行算法研究,使用C#进行工程实现比较合适,目前出现这样的情况,有一个数组,经过某种超分辨算法得到的数据点很稀疏,而且峰区域变得又高又细的。所以需要对该区域求和,就涉及到了峰位的确定,进而进行峰区域的确定,这里要注意,必须先确定峰位,再谷位,进而峰区域。
matlab实现
matlab实现算法的思路为
1、基于局部极值算法从原始数据数组获取局部极值数组(极大值,极小值,极大值索引,极小值索引);
2、极大值降序排列;
代码为
clc clear xin=load('xxx.mat'); [xmax,imax,xmin,imin] = TIAOXUAN( xin );
函数TIAOXUAN为
function [xmax,imax,xmin,imin] = TIAOXUAN( y_after_gold ) M=y_after_gold; jmin=0; jmax=0; for j=2:length(M)-1 if (M(j)>M(j-1))&&(M(j)>M(j+1)) jmax=jmax+1; max_index(jmax)=j; max_value(jmax)=M(j); end; if (M(j)<M(j-1))&&(M(j)<M(j+1)) jmin=jmin+1; min_index(jmin)=j; min_value(jmin)=M(j); end; end; xmax=max_value; xmin=min_value; imax=max_index; imin=min_index; [temp,inmax] = sort(-xmax); clear temp xmax = xmax(inmax); imax = imax(inmax); [xmin,inmin] = sort(xmin); imin = imin(inmin); end
3、查找与极大值的索引最相邻的两个极小值索引,确定峰区域;
4、峰区域中原始数据数组求和;
代码为
geshu=5; x1=zeros(geshu,2); y1=zeros(geshu,1); imax_new=imax(1:geshu); for i=1:geshu x1(i,:)=zhidao_nearest(imin,imax_new(i)); x1(i,:)=sort(x1(i,:),2); y1(i,:)=sum(xin(x1(i,1):x1(i,2))); end
函数zhidao_nearest为
function y=zhidao_nearest(A,b) [Asort,index]=sort(abs(A-b)); y(:,1)=A(index(1)); y(:,2)=A(index(2)); end
画图显示
output=[imax_new',x1,y1]; for i=1:geshu plot(output(i,2):output(i,3),xin(output(i,2):output(i,3))) hold on str=[repmat(' 峰位:',1,1) num2str(output(i,1)) char(13,10)' repmat(',峰高:',1,1) num2str(xin(output(i,1))) char(13,10)' repmat(',峰面积:',1,1) num2str(output(i,4)) char(13,10)' repmat(',峰区域:',1,1) strcat(num2str(output(i,2)),':',num2str(output(i,3)))]; text(output(i,1),xin(output(i,1)),cellstr(str)) end hold off
结果为
C#改写
C#改写存在比较多的难题,但是可以慢慢解决,下面一步一步开讲(待续)
涉及到的内容有
1、c#二维数组排序
2、。。。
直接上代码
using MathNet.Numerics.LinearAlgebra; using System; using System.Data; using System.IO; using System.Linq; namespace ConsoleApplication1 { internal class Program { /// <summary> /// 峰面积计算函数 /// </summary> /// <param name="value_index"></param> /// <param name="content"></param> /// <param name="y_after_gold"></param> /// <param name="geshu"></param> /// <param name="start"></param> /// <param name="interval"></param> /// <returns></returns> private static double[,] area_calc(double[,] value_index, Info content, double[] y_after_gold, int geshu, int start, int interval) { int[] xvalue_min = new int[geshu];//qujian left int[] xvalue_max = new int[geshu];//qujian right double[] yvalue = new double[geshu];//yvalue is peakarea double[] imax_new = new double[geshu]; //double[] value_index_min=new double[value_index.GetLength(0)]; //Generate.Map(value_index_min, x => x + 1.0); for (int i = 0; i < geshu; i++) { imax_new[i] = value_index[i, 0]; zhaodao_nearest(content.min_index, imax_new[i], out xvalue_min[i], out xvalue_max[i]); //int newlength = Convert.ToInt32(xvalue_max[i] - xvalue_min[i] + 1); //double[] temp = new double[newlength]; //foreach(var test in temp) //for (int j = xvalue_min[i]; j < xvalue_max[i]+1; j++) { for (int j = (xvalue_min[i] - start) / interval; j < (xvalue_max[i] - start) / interval + 1; j++) { //Console.WriteLine(yvalue[i]); yvalue[i] += y_after_gold[j]; //Console.WriteLine(xvalue_min[i]); } } double[,] result = new double[yvalue.Length, 4]; for (int i = 0; i < result.GetLength(0); i++) { result[i, 0] = imax_new[i]; result[i, 1] = yvalue[i]; result[i, 2] = xvalue_min[i]; result[i, 3] = xvalue_max[i]; } return result; } /// <summary> /// 主函数 /// </summary> /// <param name="args"></param> private static void Main(string[] args) { int start = 215; int interval = 15; double[] y_after_gold = loadspec("解谱结果"); //int start = 1;201 //int interval = 1; //double[] y_after_gold = loadspec("spectest"); Info content = TIAOXUAN(y_after_gold, start, interval); //string specname = "test.txt"; //try { File.Delete(@specname); } //catch { } double[,] value_index = new double[y_after_gold.Length, 4]; for (int i = 0; i < value_index.GetLength(0); i++) { value_index[i, 0] = content.max_index[i]; value_index[i, 1] = content.max_value[i]; value_index[i, 2] = content.min_index[i]; value_index[i, 3] = content.min_value[i]; } Sort(value_index, 1, "DESC"); //print2screen(value_index); savespec(value_index, "peak_location.txt"); //TODO:还未实现按列输出数据 double[,] result = area_calc(value_index, content, y_after_gold, 5, start, interval); savespec(result, "peak_area.txt"); //print2screen(result); } /// <summary> /// 寻找峰谷 /// </summary> /// <param name="arr"></param> /// <param name="number"></param> /// <param name="min"></param> /// <param name="max"></param> public static void zhaodao_nearest(double[] arr, double number, out int min, out int max) { var list = arr.ToList(); //将数组转换成List<T> list.Sort(); //排序 var l1 = list.Where(i => i < number).ToList(); //获取所有比检索值小的值 min = Convert.ToInt32(l1.Count == 0 ? 0 : l1.Max()); var l2 = list.Where(i => i > number).ToList();//获取所有比检索值大的值 max = Convert.ToInt32(l2.Count == 0 ? 0 : l2.Min()); } private static void print2screen(double[,] values) { for (int i = 0; i < values.GetLength(0); i++) { for (int k = 0; k < values.GetLength(1); k++) { Console.Write(values[i, k]); Console.Write("\t"); } Console.Write("\n"); } Console.WriteLine("\n"); } /// <summary> /// 用类定义数据类型 /// </summary> public class Info { public double[] max_index; public double[] min_index; public double[] min_value; public double[] max_value; } /// <summary> /// 冒泡排序,用于一维数组,这里没用到 /// </summary> /// <param name="a"></param> public static void MPPX(double[] a) { for (int i = 0; i < a.Length - 1; i++) { for (int j = 0; j < a.Length - 1 - i; j++) { if (a[j] < a[j + 1]) { double temp = a[j]; a[j] = a[j + 1]; a[j + 1] = temp; } } } } /// <summary> /// 二维数组排序 /// </summary> /// <typeparam name="T"></typeparam> /// <param name="array"></param> /// <param name="sortCol"></param> /// <param name="order"></param> private static void Sort<T>(T[,] array, int sortCol, string order) { int colCount = array.GetLength(1), rowCount = array.GetLength(0); if (sortCol >= colCount || sortCol < 0) throw new System.ArgumentOutOfRangeException("sortCol", "The column to sort on must be contained within the array bounds."); DataTable dt = new DataTable(); // Name the columns with the second dimension index values, e.g., "0", "1", etc. for (int col = 0; col < colCount; col++) { DataColumn dc = new DataColumn(col.ToString(), typeof(T)); dt.Columns.Add(dc); } // Load data into the data table: for (int rowindex = 0; rowindex < rowCount; rowindex++) { DataRow rowData = dt.NewRow(); for (int col = 0; col < colCount; col++) rowData[col] = array[rowindex, col]; dt.Rows.Add(rowData); } // Sort by using the column index = name + an optional order: DataRow[] rows = dt.Select("", sortCol.ToString() + " " + order); for (int row = 0; row <= rows.GetUpperBound(0); row++) { DataRow dr = rows[row]; for (int col = 0; col < colCount; col++) { array[row, col] = (T)dr[col]; } } dt.Dispose(); } /// <summary> /// 寻找局部极值,包括极大、极小、极大的索引、极小的索引 /// </summary> /// <param name="M"></param> /// <param name="start"></param> /// <param name="interval"></param> /// <returns></returns> private static Info TIAOXUAN(double[] M, int start, int interval) { Info info = new Info(); int jmin = 0, jmax = 0; info.max_index = new double[M.Length]; info.min_index = new double[M.Length]; info.min_value = new double[M.Length]; info.max_value = new double[M.Length]; for (int j = 1; j < M.Length - 2; j++) { jmax++; if ((M[j] > M[j - 1]) && (M[j] > M[j + 1])) { info.max_index[jmax] = j * interval + start; info.max_value[jmax] = M[j]; } else { info.max_index[jmax] = 0; } jmin++; if ((M[j] < M[j - 1]) && (M[j] < M[j + 1])) { info.min_index[jmin] = j * interval + start; info.min_value[jmax] = M[j]; } else { info.min_index[jmin] = 0; } } return info; } /// <summary> /// 载入数据 /// </summary> /// <param name="specname"></param> /// <returns></returns> private static double[] loadspec(string specname) { var M = Matrix<double>.Build; var V = Vector<double>.Build; string[] content = File.ReadAllLines(@specname); double[] yvalue = new double[content.Length]; double[] axesx = new double[content.Length]; for (int i = 1; i < content.Length; i++) { //axesx[i] = i; yvalue[i] = Convert.ToDouble(content[i]); } return yvalue; } /// <summary> /// 保存二维数组数据 /// </summary> /// <param name="content"></param> /// <param name="specname"></param> private static void savespec(double[,] content, string specname) { FileStream fs = new FileStream(@specname, FileMode.Append); StreamWriter sw = new StreamWriter(fs); //method one //for (int i = 1; i < content.Length; i++) //{ // string str = content[i].ToString() + ","; // str.TrimEnd(); // //sw.Write(content[i] + " ");//Convert.ToString(hist) // sw.Write(str.TrimEnd(','));//Convert.ToString(hist) //} //method two //string str = "["; //foreach (double test in content) // str += test.ToString() + ","; // //str += string.Format("{0:E2}",test)+","; //str = str.TrimEnd(',') + "];"; //sw.Write(str);//Convert.ToString(hist) //sw.Write("\n"); //method three for (int i = 0; i < content.GetLength(0); i++) { for (int k = 0; k < content.GetLength(1); k++) { sw.Write(content[i, k]); sw.Write("\t"); } sw.Write("\n"); } sw.WriteLine("\n"); sw.Close(); sw.Dispose(); sw = null; } } }
实线的效果
参考文献:
1、脚本之家:C#实现对二维数组排序的方法
2、ITPUB网站.NET技术的博客:C# 实现二维数组的排序算法(代码)
3、Lic.的matlab局部极值算法(代码)