how to implement savitzky golay filter without using inbuilt functions
13 views (last 30 days)
Show older comments
code for savitzky golay filter without using sgolayfilt() to perform smoothing and detect peaks in a signal
0 Comments
Answers (5)
John D'Errico
on 14 Apr 2017
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 Comment
Vrushabh Bhangod
on 21 May 2018
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
on 14 Apr 2017
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 Comments
Image Analyst
on 22 May 2018
Well, it's not necessarily the best. If the curve goes upward or bends around, you might use middleX = mean(x).
Vrushabh Bhangod
on 24 May 2018
Edited: Vrushabh Bhangod
on 24 May 2018
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]);
0 Comments
Zeus
on 28 May 2018
Edited: Zeus
on 28 May 2018
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);
0 Comments
Shahzad
on 24 Sep 2022
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);
1 Comment
Image Analyst
on 24 Sep 2022
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!