Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.DataImporter/HeuristicLab.DataImporter.Command/ChangeValues/FilterSavitzkyGolayCommand.cs @ 6134

Last change on this file since 6134 was 6134, checked in by gkronber, 13 years ago

#1471: added plugin for DbExplorer interfaces, deleted .resx files, set svn:ignore properties, and added license header

File size: 6.5 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2011 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Linq;
25using System.Text;
26using System.Xml;
27using HeuristicLab.DataImporter.Data;
28using HeuristicLab.DataImporter.Data.CommandBase;
29using HeuristicLab.DataImporter.Data.Model;
30using HeuristicLab.DataImporter.Command.View;
31using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
32
33namespace HeuristicLab.DataImporter.Command {
34  [StorableClass]
35  [ViewableCommandInfoAttribute("Savitzky Golay", 1, ColumnGroupState.DoubleColumnSelected, "Change Values", Position = 2,
36    OptionsView = typeof(FilterSavitzkyGolayCommandView))]
37  public class FilterSavitzkyGolayCommand : ColumnGroupCommandWithAffectedColumnsBase {
38    private Dictionary<int, DoubleColumn> oldColumns;
39    private ICollection<int> oldSortedColumnIndices;
40
41    public int WindowLeft { get; set; }
42
43    public int WindowRight { get; set; }
44
45    public int Order { get; set; }
46
47    private FilterSavitzkyGolayCommand()
48      : base(null, string.Empty, null) {
49      oldColumns = new Dictionary<int, DoubleColumn>();
50    }
51
52    public FilterSavitzkyGolayCommand(DataSet dataSet, string columnGroupName, int[] affectedColumns)
53      : base(dataSet, columnGroupName, affectedColumns) {
54      oldColumns = new Dictionary<int, DoubleColumn>();
55      this.WindowLeft = -16;
56      this.WindowRight = 16;
57      this.Order = 2;
58    }
59
60    public override string Description {
61      get { return "Filter Column - Savitzky Golay"; }
62    }
63
64    public override void Execute() {
65      base.Execute();
66      if (this.WindowLeft >= WindowRight)
67        throw new CommandExecutionException("Window size for filter commands must no be smaller than 1.", this);
68      foreach (int col in AffectedColumns) {
69        DoubleColumn doubleCol = ColumnGroup.GetColumn(col) as DoubleColumn;
70        if (doubleCol == null) throw new CommandExecutionException("Filtering is only supported for double columns.", this);
71        if (doubleCol.ContainsNullValues) throw new CommandExecutionException("Filtering is not supported for double columns with missing values.", this);
72      }
73      DoubleColumn column;
74      oldSortedColumnIndices = new List<int>(ColumnGroup.SortedColumnIndexes);
75      foreach (int col in AffectedColumns) {
76        if (ColumnGroup.GetColumn(col) is DoubleColumn) {
77          column = (DoubleColumn)ColumnGroup.GetColumn(col);
78          oldColumns.Add(col, column);
79          ColumnGroup.ReplaceColumn(col, CalcNewColumn(column, this.WindowLeft, this.WindowRight, this.Order));
80        }
81      }
82      ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
83      ColumnGroup.FireChanged();
84      ColumnGroup = null;
85    }
86
87    public override void UndoExecute() {
88      base.UndoExecute();
89      foreach (KeyValuePair<int, DoubleColumn> pair in oldColumns)
90        ColumnGroup.ReplaceColumn(pair.Key, pair.Value);
91
92      ColumnGroup.SortedColumnIndexes = oldSortedColumnIndices;
93      oldSortedColumnIndices = null;
94      oldColumns.Clear();
95      ColumnGroup.FireChanged();
96      ColumnGroup = null;
97    }
98
99    private DoubleColumn CalcNewColumn(DoubleColumn oldColumn, int windowLeft, int windowRight, int order) {
100      double[] c;
101      SavitzkyGolay(Math.Abs(windowLeft), Math.Abs(windowRight), order, out c);
102
103      DoubleColumn newCol = (DoubleColumn)oldColumn.CreateCopyOfColumnWithoutValues();
104      double[] data = oldColumn.ValuesEnumerable.Cast<double>().ToArray();
105      double[] ans = new double[oldColumn.TotalValuesCount];
106      alglib.convr1d(data, data.Length, c, c.Length, out ans);
107      for (int i = Math.Abs(windowRight); i < ans.Length - Math.Abs(windowLeft); i++) {
108        newCol.AddValue(ans[i]);
109      }
110      newCol.SortOrder = oldColumn.SortOrder;
111      return newCol;
112    }
113
114    /// <summary>
115    /// Calculates coefficients for Savitzky-Golay filtering. (Numerical Recipes, page 769), one important change is that the coefficients are returned in normal order instead of wraparound order
116    /// </summary>
117    /// <param name="nl">number of samples to the left</param>
118    /// <param name="nr">number of samples to the right</param>
119    /// <param name="m">order of the polynomial to fit</param>
120    /// <param name="c">resulting coefficients for convolution, in correct order (t-nl, ... t-1, t+0, t+1, ... t+nr)</param>
121    private void SavitzkyGolay(int nl, int nr, int order, out double[] c) {
122      int np = nl + nr + 1;
123      int ld = 0;
124      int m = order;
125
126      int j, k, imj, ipj, kk, mm;
127      double fac = 0;
128      double sum = 0;
129      if (np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m) throw new ArgumentException();
130
131      int[] indx = new int[m + 1];
132      double[,] a = new double[m + 1, m + 1];
133      double[] b = new double[m + 1];
134      c = new double[np];
135
136      for (ipj = 0; ipj <= (m << 1); ipj++) {
137        sum = (ipj > 0 ? 0.0 : 1.0);
138        for (k = 1; k <= nr; k++) sum += Math.Pow((double)k, (double)ipj);
139        for (k = 1; k <= nl; k++) sum += Math.Pow((double)-k, (double)ipj);
140        mm = Math.Min(ipj, 2 * m - ipj);
141        for (imj = -mm; imj <= mm; imj += 2)
142          a[(ipj + imj) / 2, (ipj - imj) / 2] = sum;
143      }
144      for (j = 0; j < m + 1; j++) b[j] = 0;
145      b[ld] = 1.0;
146      alglib.densesolverreport rep;
147      int info;
148      double[] x = new double[b.Length];
149      alglib.rmatrixsolve(a, b.Length, b, out info, out rep, out x);
150
151      for (kk = 0; kk < np; kk++) c[kk] = 0.0;
152      for (k = -nl; k <= nr; k++) {
153        sum = x[0];
154        fac = 1.0;
155        for (mm = 1; mm <= m; mm++) sum += x[mm] * (fac *= k);
156        kk = k + nl;
157        c[kk] = sum;
158      }
159    }
160  }
161}
Note: See TracBrowser for help on using the repository browser.