how to implement savitzky golay filter without using inbuilt functions
13 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
code for savitzky golay filter without using sgolayfilt() to perform smoothing and detect peaks in a signal
0 commentaires
Réponses (5)
John D'Errico
le 14 Avr 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 commentaire
Vrushabh Bhangod
le 21 Mai 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
le 14 Avr 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 commentaires
Image Analyst
le 22 Mai 2018
Well, it's not necessarily the best. If the curve goes upward or bends around, you might use middleX = mean(x).
Vrushabh Bhangod
le 24 Mai 2018
Modifié(e) : Vrushabh Bhangod
le 24 Mai 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 commentaires
Zeus
le 28 Mai 2018
Modifié(e) : Zeus
le 28 Mai 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 commentaires
Shahzad
le 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 commentaire
Image Analyst
le 24 Sep 2022
Voir également
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!