Как работает функция sgolay в Matlab R2013a?

У меня есть вопрос о функции sgolay в Matlab R2013a. В моей базе данных 165 спектров с 2884 переменными, и я хотел бы взять из них первую и вторую производные. Как я могу определить входы K и F до sgolay?

Ниже приведен пример:

sgolay используется для сглаживания зашумленной синусоиды и сравнения полученных первой и второй производных с первой и второй производными, вычисленными с помощью diff. Обратите внимание, как использование diff усиливает шум и приводит к бесполезным результатам.

K = 4;                 % Order of polynomial fit
F = 21;                % Window length
[b,g] = sgolay(K,F);   % Calculate S-G coefficients

dx = .2;
xLim = 200;
x = 0:dx:xLim-1;

y = 5*sin(0.4*pi*x)+randn(size(x));  % Sinusoid with noise

HalfWin  = ((F+1)/2) -1;
for n = (F+1)/2:996-(F+1)/2,
  % Zero-th derivative (smoothing only)
  SG0(n) =   dot(g(:,1), y(n - HalfWin: n + HalfWin));

  % 1st differential
  SG1(n) =   dot(g(:,2), y(n - HalfWin: n + HalfWin));

  % 2nd differential
  SG2(n) = 2*dot(g(:,3)', y(n - HalfWin: n + HalfWin))';
end

SG1 = SG1/dx;         % Turn differential into derivative
SG2 = SG2/(dx*dx);    % and into 2nd derivative

% Scale the "diff" results
DiffD1 = (diff(y(1:length(SG0)+1)))/ dx;    
DiffD2 = (diff(diff(y(1:length(SG0)+2)))) / (dx*dx);

subplot(3,1,1);
plot([y(1:length(SG0))', SG0'])
legend('Noisy Sinusoid','S-G Smoothed sinusoid')

subplot(3, 1, 2);
plot([DiffD1',SG1'])
legend('Diff-generated 1st-derivative', 'S-G Smoothed 1st-derivative')

subplot(3, 1, 3);
plot([DiffD2',SG2'])
legend('Diff-generated 2nd-derivative', 'S-G Smoothed 2nd-derivative')

matlab filtering smoothing derivative
person user3689163    schedule 29.05.2014    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Взятие производных в изначально зашумленном процессе. Таким образом, если у вас уже есть некоторый шум в ваших данных, он действительно будет увеличиваться по мере того, как вы берете производные более высокого порядка. Savitzky-Golay — очень полезный способ объединения сглаживания и дифференцирования в одну операцию. Это общий метод, и он вычисляет производные в произвольном порядке. Однако есть компромиссы. Существуют и другие специальные методы для данных с определенной структурой.

Что касается вашего заявления, у меня нет конкретных ответов. Многое зависит от характера данных (частота дискретизации, коэффициент шума и т. д.). Если вы используете слишком много сглаживания, вы размажете свои данные или создадите алиасинг. То же самое, если вы подгоняете данные, используя коэффициенты полинома высокого порядка, K. В вашем демонстрационном коде вы также должны построить аналитические производные функции sin. Затем поиграйте с различным количеством входного шума и сглаживающими фильтрами. Такой инструмент с известными точными ответами может быть полезен, если вы можете аппроксимировать аспекты ваших реальных данных. На практике я стараюсь использовать как можно меньше сглаживания, чтобы производить не слишком шумные производные. Часто это означает полином третьего порядка (K = 3) и размер окна F как можно меньше.

Так что да, многие предлагают вам смотреть на глаз для настройки этих параметров. Однако совсем недавно было проведено несколько исследований по автоматическому выбору коэффициентов: оптимальных фильтров Савицкого-Голея (2013). Существуют также альтернативы методу Савицкого-Голея, например, эта статья на основе регуляризация, но вам может понадобиться реализовать их самостоятельно в Matlab.

Кстати, некоторое время назад я написал небольшую замену для sgolay. Как и вам, мне нужен был только второй вывод, фильтры дифференцирования, G, так что это все, что он вычисляет. Эта функция также быстрее (примерно в 2–4 раза):

function G=sgolayfilt(k,f)
%SGOLAYFILT  Savitzky-Golay differentiation filters
s = vander(0.5*(1-f):0.5*(f-1));
S = s(:,f:-1:f-k);
[~,R] = qr(S,0);
G = S/R/R';

Полная версия этой функции с проверкой ввода доступна на моем GitHub.

person horchler    schedule 30.05.2014
comment
Спасибо за пояснение и указание статьи hockler, однако все же есть сомнения в части вывода, который вы предложили. Я не близок к этой части программирования, вы могли бы объяснить мне, что это значит и как каждая формула получает мои 165 предварительно обработанных спектров (первая и вторая производная). Извините за невежество, а также за то, что я пишу, не могу писать по английскому праву. - person user3689163; 02.06.2014
comment
Спасибо horchler за помощь, проблема решена. Мне нужно было немного больше изучить sgolay. Извините что-нибудь. - person user3689163; 04.06.2014