arntanguy / gram_savitzky_golay

C++ Implementation of Savitzky-Golay filtering based on Gram polynomials
BSD 2-Clause "Simplified" License
100 stars 29 forks source link

Simple smoothing? #11

Open yrzhanov opened 2 years ago

yrzhanov commented 2 years ago

Hello, I think that SG filter may perform better in my case than the simple moving average. I have 1D static equispaced data and I just need it smoothed - with a filter length found by trial and error. Could you please provide an example that would return smoothed Y value for any required X? I expected something like: filter = MakeFilter(length) filter.SmoothData(rawData = rawX and rawY) smoothY = filter.SmoothValue(rawX)

but the filter takes all the raw data as an agrument.

Thank you, Yuri

heisenbergpxh commented 6 months ago

It might be little late but if you want to generate the same output as for instance scipy's savgol_filter

from scipy.signal import savgol_filter

x = [.1, .2, .3, .4, .5, .6, .7]

from scipy.signal import savgol_filter

savgol_filter(x, 5, 2, deriv=1, mode='constant')

which returns [ 0.08 0.1 0.1 0.1 0.1 -0.06 -0.16]. With this C++ implementation you have to perform a moving window over your data which have to be extended at the edges with 0.


const int m = 2; // Window size is 2*m+1
const int n = 2; // Polynomial Order
const int t = 0; // Evaluate the polynomial in the center (index 2 of the 5 entry long vector) 
const int d = 1; // First order derivation

gram_sg::SavitzkyGolayFilter filter(m, t, n, d);

std::vector<double> filtered_data;
int half_window_size = (2 * m + 1) / 2;

std::vector<double> data = {.1, .2, .3, .4, .5, .6, .7};

data.insert(data.begin(), half_window_size, 0);
data.insert(data.end(), half_window_size, 0);

// data = {0, 0, .1, .2, .3, .4, .5, .6, .7, 0, 0};

for (int i = 0; i < data.size(); ++i) {
  std::vector<double> window;

  int lower_idx = std::max(0, i - half_window_size);
  int upper_idx = std::min(i + half_window_size, static_cast<int>(data.size()) - 1);

  for (int j = lower_idx; j <= upper_idx; ++j) {
      window.push_back(data[j]);
  }

  filtered_data.push_back(filter.filter(window));
}

filtered_data.erase(filtered_data.begin(), filtered_data.begin() + half_window_size);
filtered_data.erase(filtered_data.end() - half_window_size, filtered_data.end());

for (const auto& elem : filtered_data) {
  std::cout << elem << " ";
}

Returns exactly the same 0.08 0.1 0.1 0.1 0.1 -0.06 -0.16

JimmyDaSilva commented 5 months ago

Thanks a lot @heisenbergpxh . I am also trying to get the same results as Scipy.

Do you have an equivalent for Scipy's interp mode ?

yrzhanov commented 5 months ago

Thank you, good algorithm is never late. And you never know when you need it.