how to implement savitzky golay filter without using inbuilt functions

code for savitzky golay filter without using sgolayfilt() to perform smoothing and detect peaks in a signal

Answer 1

[John D'Errico]

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.

[Vrushabh Bhangod]

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.

Answer 2

[Image Analyst]

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.

[Vrushabh Bhangod]

Sir/Madam, Can you kindly explain the evaluate the polynomial statement...I have constructed the code until here and I am unable to understand what the further statement means. Please reply

[Image Analyst]

for k = 1 : length(x)-windowWidth-1
  x = k : (k + windowWidth + 1) % Extract x within the window.  
  % Do the fit with x and y that are just in the window only.
  coefficients = polyfit(x, y, 2)
  % Get the fitted value at the middle of the window.
  middleX = (x(1) + x(end))/2;
  middleValue = polyval(coefficients, middleX);

and so on...

[Vrushabh Bhangod]

Thanks a lot,i had evaluated middlex = x(mid) :(

[Image Analyst]

Well, it's not necessarily the best. If the curve goes upward or bends around, you might use middleX = mean(x).

Answer 3

[Vrushabh Bhangod]

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]);

Answer 4

[Zeus]

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);

Answer 5

[Shahzad]

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);

[Image Analyst]

@Shahzad, didn't you just repost @Zeus's code? Why? Do you want me to delete it for you?