做数据分析算法,使用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局部极值算法(代码)