数组数据分析算法中峰区域的确定

做数据分析算法，使用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("	");
                }
                Console.Write("
");
            }
            Console.WriteLine("
");
        }

        /// <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("
");

            //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("	");
                }
                sw.Write("
");
            }
            sw.WriteLine("
");

            sw.Close();
            sw.Dispose();
            sw = null;
        }
    }
}

　　实线的效果

 

参考文献：

1、脚本之家：C#实现对二维数组排序的方法

2、ITPUB网站.NET技术的博客：C# 实现二维数组的排序算法(代码) 

3、Lic.的matlab局部极值算法（代码）