how to implement savitzky golay filter without using inbuilt functions
조회 수: 13 (최근 30일)
이전 댓글 표시
code for savitzky golay filter without using sgolayfilt() to perform smoothing and detect peaks in a signal
댓글 수: 0
답변 (5개)
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
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
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
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
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]);
댓글 수: 0
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);
댓글 수: 0
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);
댓글 수: 1
Image Analyst
2022년 9월 24일
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!