using FlyADBase; using GalaSoft.MvvmLight.Command; using MathNet.Numerics.LinearAlgebra.Double; using System; using System.Collections.Generic; using System.Collections.ObjectModel; using System.ComponentModel; using System.Linq; using System.Text; using System.Threading.Tasks; namespace Flyad7_WPF { public class WdConvolutionVm : INotifyPropertyChanged { public event PropertyChangedEventHandler PropertyChanged; public RelayCommand CreateHolderCmd { get; private set; } public RelayCommand CalConvCmd { get; private set; } public RelayCommand CalConv2Cmd { get; private set; } public RelayCommand UpdateConvCmd { get; private set; } /// /// 小托辊不锈钢支架 宽度单位 mm, 默认15mm /// public int HolderWidth { get; set; } = 15; /// /// 编码器1 的比例 mm/脉冲 /// public double Mmpp { get; set; } = 0.16; /// /// 小托辊不锈钢支架 中间grid位置 /// public int HolderMidGrid { get; set; } = 40; /// /// 小托辊不锈钢支架 开始grid位置 /// public int HolderBeginGrid { get; set; } = 20; /// /// 小托辊不锈钢支架 结束grid位置 /// public int HolderEndGrid { get; set; } = 60; public int HolderLenGrid => HolderEndGrid - HolderBeginGrid + 1; /// /// 空气AD /// public int AirAd { get; set; } = 25900; /// /// 小托辊不锈钢支架 grid图 /// public int[] HolderDatas { get; set; } /// /// 小托辊不锈钢支架 滤波后的 grid图, 就是AD盒采集的 /// public int[] FilterDatas { get; set; } /// /// 把滤波后的数据,反卷积后的 grid图 /// public int[] DeconvDatas { get; set; } /// /// 卷积核 大小,单位grid /// public ObservableCollection ConvCore { get; } = new ObservableCollection(); /// /// 卷积核 大小,单位grid /// public int CoreLen { get; set; } = 17; /// /// DeconvDatas 与 HolderDatas 的相关性 /// public double R { get; set; } /// /// FilterDatas 与 HolderDatas 的相关性 /// public double HolderR { get; set; } System.Windows.Forms.DataVisualization.Charting.Chart chart0; System.Windows.Forms.DataVisualization.Charting.Chart chart1; System.Windows.Forms.DataVisualization.Charting.Chart chart2; FlyAD7 flyad; public WdConvolutionVm() { CreateHolderCmd = new RelayCommand(CreateHolder); CalConvCmd = new RelayCommand(CalConv); CalConv2Cmd = new RelayCommand(CalConv2); UpdateConvCmd = new RelayCommand(UpdateConv); } public void Init( FlyAD7 flyad, System.Windows.Forms.DataVisualization.Charting.Chart chart0, System.Windows.Forms.DataVisualization.Charting.Chart chart1, System.Windows.Forms.DataVisualization.Charting.Chart chart2, int[] filterDatas) { this.flyad = flyad; this.chart0 = chart0; this.chart1 = chart1; this.chart2 = chart2; FilterDatas = filterDatas; DrawFilter(); } private async void DrawFilter() { await Draw0(0, FilterDatas); } async Task Draw0(int index, int[] data) { var chart = chart0; await Task.Yield(); //画图 await Task.Factory.StartNew(() => { App.Current.Dispatcher.Invoke(() => { System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index]; series.Points.Clear(); for (int i = 0; i < data.Count(); i++) { if (Misc.MyBase.ISVALIDATA(data[i])) series.Points.AddXY(i, data[i]); } //if (series.Points.Count() > 0) //{ // double max = series.Points.Max(p => p.YValues[0]); // double min = series.Points.Min(p => p.YValues[0]); // if (max != min) // { // double range = (max - min) / 10; // chart.ChartAreas[0].AxisY.Minimum = min - range; // chart.ChartAreas[0].AxisY.Maximum = max + range; // } //} }); }); } async Task Draw2(int index, double r) { var chart = chart2; //画图 await Task.Yield(); await Task.Factory.StartNew(() => { App.Current.Dispatcher.Invoke(() => { System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index]; series.Points.AddY(r * 10000); double max = series.Points.Max(p => p.YValues[0]); double min = series.Points.Min(p => p.YValues[0]); if (max != min) { double range = (max - min) / 10; chart.ChartAreas[0].AxisY.Minimum = min - range; chart.ChartAreas[0].AxisY.Maximum = max + range; } }); }); } async Task Draw1(int index,IEnumerable data) { var chart = chart1; //画图 await Task.Yield(); await Task.Factory.StartNew(() => { App.Current.Dispatcher.Invoke(() => { System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index]; series.Points.Clear(); for (int i = 0; i < data.Count(); i++) { series.Points.AddXY(i, data.ElementAt(i)); } double max = series.Points.Max(p => p.YValues[0]); double min = series.Points.Min(p => p.YValues[0]); if (max != min) { double range = (max - min) / 10; chart.ChartAreas[0].AxisY.Minimum = min - range; chart.ChartAreas[0].AxisY.Maximum = max + range; } }); }); } private async void CalConv() { //解卷积核 // Holder * M = Filter; // A * X = B; DenseMatrix matrixA = new DenseMatrix(HolderLenGrid, CoreLen); DenseVector vectorB = new DenseVector(HolderLenGrid); for (int i = 0; i < HolderLenGrid; i++) { DenseVector vector = new DenseVector(CoreLen); for (int j = 0; j < CoreLen; j++) { vector[j] = HolderDatas[HolderBeginGrid + i + j - CoreLen / 2]; } matrixA.SetRow(i, vector); vectorB[i] = FilterDatas[HolderBeginGrid + i]; } //A~AX = A~B DenseMatrix ATA = (DenseMatrix)matrixA.TransposeThisAndMultiply(matrixA); DenseVector ATB = (DenseVector)matrixA.TransposeThisAndMultiply(vectorB); //DenseVector X = (DenseVector)matrixA.LU().Solve(vectorB); DenseVector X = (DenseVector)ATA.LU().Solve(ATB); ConvCore.Clear(); //导出到 ConvCore for (int i = 0; i < CoreLen; i++) { ConvCore.Add(new CoreValue() { Value = X[i] }); } await Draw1(0, ConvCore.Select(c=>c.Value)); var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value)); await Draw0(2, filterdatas2); R = Misc.MyMath.Correl(FilterDatas, filterdatas2); Draw2Clear(); await Draw2(0, R); } async void Draw2Clear() { var chart = chart2; await Task.Yield(); //画图 App.Current.Dispatcher.Invoke(() => { System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[0]; series.Points.Clear(); //chart.ChartAreas[0].AxisX.Minimum = double.NaN; //chart.ChartAreas[0].AxisX.Maximum = double.NaN; }); } private async void CalConv2() { double step = 0.00001; int decimals = 5; //在 ConvCore 基础上微调 //必须为单数 if (ConvCore.Count % 2 != 1) ConvCore.RemoveAt(ConvCore.Count - 1); //必须都大于等于0 //ConvCore 必须是对称的 for (int i = 0; i < ConvCore.Count() / 2; i++) { if (ConvCore[i].Value < 0) ConvCore[i].Value = 0; ConvCore[i].Value = Math.Round(ConvCore[i].Value, decimals); ConvCore[ConvCore.Count - 1 - i].Value = ConvCore[i].Value; } while (true) { var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value)); double r0 = Misc.MyMath.Correl(FilterDatas, filterdatas2); bool hasChanged = false; double[] ks = ConvCore.Select(c => c.Value).ToArray(); double r = r0; await Task.Factory.StartNew(() => { //每个参数改变 step, 确定方向 for (int i = 0; i <(ks.Count()/2+1); i++) { double curr_step = step; double r1; double r2; while (true) { ks[i] = ConvCore[i].Value + curr_step; ks[ks.Count() - 1 - i] = ks[i]; filterdatas2 = Conv(HolderDatas, ks); r1 = Misc.MyMath.Correl(FilterDatas, filterdatas2); ks[i] = ConvCore[i].Value - curr_step; ks[ks.Count() - 1 - i] = ks[i]; filterdatas2 = Conv(HolderDatas, ks); r2 = Misc.MyMath.Correl(FilterDatas, filterdatas2); if ((r1 == r0) || (r2 == r0)) { //居然没有变化,把step翻倍 curr_step *= 1.5; } else { //有变化了 break; } } if (r1 > r0) { //方向正确 ks[i] = ConvCore[i].Value + curr_step; ks[ks.Count() - 1 - i] = ks[i]; hasChanged = true; r = r1; continue; } else if (r2 > r0) { //方向正确 ks[i] = ConvCore[i].Value - curr_step; ks[ks.Count() - 1 - i] = ks[i]; hasChanged = true; r = r2; continue; } else { //方向错误,回到以前 ks[i] = ConvCore[i].Value; } } }); R = r; await Draw2(0, r); if (!hasChanged) { //已经计算完毕 break; } else { for (int i = 0; i < ks.Count(); i++) ConvCore[i].Value = ks[i]; await Draw1(0, ConvCore.Select(c => c.Value)); } } UpdateConv(); } private async void UpdateConv() { await Draw1(0, ConvCore.Select(c => c.Value)); var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value)); await Draw0(2, filterdatas2); R = Misc.MyMath.Correl(FilterDatas, filterdatas2); await Draw2(0, R); } int[] Conv(int[] a, IEnumerable core) { int[] b = new int[a.Length]; for (int i = 0; i < a.Length; i++) { double sum=0; double cnt=0; for (int j = 0; j < core.Count(); j++) { int index = i - core.Count() / 2 + j; if (index < 0) continue; if (index > a.Length - 1) break; if (Misc.MyBase.ISVALIDATA(a[index])){ cnt += core.ElementAt(j); sum += core.ElementAt(j) * a[index]; } } if (cnt > 0) { b[i] = (int)(sum / cnt); } else { b[i] = Misc.MyBase.NULL_VALUE; } } return b; } private async void CreateHolder() { HolderDatas = new int[FilterDatas.Count()]; //把 Target 的 测到不锈钢部分删除, 剩下求平均值 int sensorGridWidth = (int)(25 / Mmpp / flyad.PosOfGrid); int holderGridWidth = (int)(HolderWidth / Mmpp / flyad.PosOfGrid); for (int i = 0; i < HolderDatas.Count(); i++) { if (Misc.MyBase.ISVALIDATA(FilterDatas[i]))// && i>= HolderBeginGrid && i<= HolderEndGrid) HolderDatas[i] = AirAd; else HolderDatas[i] = Misc.MyBase.NULL_VALUE; } for (int i = 0; i < holderGridWidth; i++) { int index = HolderMidGrid + i - holderGridWidth / 2; HolderDatas[index] = 0; } //找到最大相关性对应的位置 HolderR = Misc.MyMath.Correl(FilterDatas, HolderDatas); await Draw0(1, HolderDatas); } } public class CoreValue :INotifyPropertyChanged { public double Value { get; set; } public event PropertyChangedEventHandler PropertyChanged; } }