how to implement savitzky golay filter without using inbuilt functions

51 次查看（过去 30 天）
preeti visweswaran
preeti visweswaran 2017-4-14
评论： Image Analyst 2022-9-24
code for savitzky golay filter without using sgolayfilt() to perform smoothing and detect peaks in a signal

请先登录，再回答此问题。

回答（5 个）

John D'Errico
John D'Errico 2017-4-14
Learn to use tools like conv or filter. They can accomplish the desired result, given the proper input. Or download a Savitsky-Golay tool from the file exchange. As I recall, there are lots of them.
  1 个评论
Vrushabh Bhangod
Vrushabh Bhangod 2018-5-21
Sir, You mean that I have to perform Discrete convolution with a fixed impulse response? I read an IEEE paper which describes SG filter that way.

请先登录，再进行评论。

Image Analyst
Image Analyst 2017-4-14
You can use a sliding window. Just march your window along your signal and extract the data in the window and fit it to a polynomial using polyfit(). Evaluate the polynomial at the center element location with polyval(), and that's your output for that location. Really pretty easy. I suggest you at least give it an attempt yourself. I think you know how to use a for loop and how to call polyfit() and polyval() so it should be trivial.
  4 个评论
显示 2更早的评论
Image Analyst
Image Analyst 2018-5-22
Well, it's not necessarily the best. If the curve goes upward or bends around, you might use middleX = mean(x).

请先登录，再进行评论。

Vrushabh Bhangod
Vrushabh Bhangod 2018-5-24
编辑：Vrushabh Bhangod 2018-5-24
if true
% code
end%%contruction of signal and addition of WG noise
n = (1:4096); % time vector
N1 = 4096;% length of signal
sig = MakeSignal('Piece-Regular',N1); %loading Piece regular Signal of length n
SNR = 10; %In dB
x = awgn(sig,SNR,'measured'); % addition of white gaussian noise
%%construction of Savitzky-Golay Filter
WinL = 15; %in samples
Ord = 3; % order of the filter
shiftL = 1; % hop size in samples
nFr = round(length(x)/shiftL); %no., of frames
WIND = zeros(WinL,nFr);
for c = 1:nFr - round(WinL/shiftL)
FB = (c-1)*shiftL+1; % beginning of the frame in samples
FE = FB + WinL -1; % ending of the frame in samples
WIND(:,c) = x(FB:FE);
end
for c = 1:nFr - round(WinL/shiftL) % computing no., of frames into windows
FB = (c-1)*shiftL+1; % beginning of the frame in samples
FE = FB + WinL -1; % ending of the frame in samples
N(:,c) = n(FB:FE);
end
adj = zeros(WinL,size(WIND,2)-size(N,2));
WIND(:,[size(N,2)+1:size(WIND,2)]) = [];
polcoeff = zeros(Ord+1,size(N,2)); % coefficients of the polynomial
polvalues = zeros(WinL,size(N,2)); % value of the function with 'p' polynomial coefficient
for c = 1:size(N,2)
t = N(:,c);
[p,s,mu] = polyfit(t,WIND(:,c),Ord);
polcoeff(:,c) = p;
polvalues(:,c) = polyval(p,t(round(WinL/2)),s,mu);
end
polvalues(2:WinL,:) = [];
polvalues = [polvalues,zeros(1,WinL)];
%%to calculate mean square error
t = sum((sig-polvalues).^2); % x is the signal with AWGN , polvalues is the recovered signal
MSE = t/N1
subplot(311)
plot(sig); ylabel('Amplitude'),xlabel('Number of samples');title('Original signal');axis([0 4096 -100 100]);
subplot(312)
plot(x);ylabel('Amplitude'),xlabel('Number of samples');title('Original signal + Noise');axis([0 4096 -100 100]);
subplot(313)
plot(polvalues);ylabel('Amplitude'),xlabel('Number of samples');title('Recovered signal');axis([0 4096 -100 100]);
Zeus
Zeus 2018-5-28
编辑：Zeus 2018-5-28
function y=sgfilter(x,ML,MR,N)
% the window size is ML+MR+1
% x is the input signal(with or without noise)
% N is the order of the polynomial that will approximate signal x in each window
% refer IEEE paper of Robert Schafer 'What is Savitzky Golay Filter?' for better understanding.
len=length(x);
xn=[zeros(1,ML),x,zeros(1,MR)];
y=zeros(1,len);
d=[-ML:MR]';
l=length(d);
A=zeros(l,N+1);
A(:,1)=1;
for i=1:N,
A(:,i+1)=d(:,1).^i;
end
H=pinv(A'*A)*A';% fliplr(H(1,:)) is actually the impulse response of the savitzky-golay filter.
for i=1:len,
in=xn(1,i:ML+MR+i);
in=in(:);
y(1,i)=H(1,:)*in;% convolution of the sgfilter's impulse response with the signal values in each window
end
figure(1);
subplot(311); plot(x);subplot(312);plot(y);subplot(313);plot(x-y);
Shahzad
Shahzad 2022-9-24
function y=sgfilter(x,ML,MR,N)
% the window size is ML+MR+1
% x is the input signal(with or without noise)
% N is the order of the polynomial that will approximate signal x in each window
% refer IEEE paper of Robert Schafer 'What is Savitzky Golay Filter?' for better understanding.
len=length(x);
xn=[zeros(1,ML),x,zeros(1,MR)];
y=zeros(1,len);
d=[-ML:MR]';
l=length(d);
A=zeros(l,N+1);
A(:,1)=1;
for i=1:N,
A(:,i+1)=d(:,1).^i;
end
H=pinv(A'*A)*A';% fliplr(H(1,:)) is actually the impulse response of the savitzky-golay filter.
for i=1:len,
in=xn(1,i:ML+MR+i);
in=in(:);
y(1,i)=H(1,:)*in;% convolution of the sgfilter's impulse response with the signal values in each window
end
figure(1);
subplot(311); plot(x);subplot(312);plot(y);subplot(313);plot(x-y);

类别

在 Help Center 和 File Exchange 中查找有关 Smoothing and Denoising 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by