Plots of Savitzky-Golay

Hi,

I have several questions regarding the Savitzky-Golay filter.

I managed to follow the example from the website here: http://uk.mathworks.com/help/signal/ref/sgolay.html

I basically want to show visually, several signals filtered on top of the noisy signal. The example given only shows for one filter and the default colours. Here are the codes that I currently used.

1) This is saved as a .mat file so I have the same noisy signal to use with different filters.

xLim = 1000; % Generate numbers from 0 to 1000.

dx = 0.1; % Increment of 0.1.

x = 0:dx:xLim-1; % Generate a series of numbers from 0 to 999 (xLim-1) with increment of 0.1.

y0 = 6*sin(0.4*pi*x); % Create sine signal.

y = awgn(y0,0.01,'measured'); % Add white Gaussian noise.

2) I load the above and then apply the filter.

N = 2; % Order of polynomial fit (2,4,6,8 required)

F = 11; % Window length (11,23,35,47,59 required)

[b,g] = sgolay(N,F); % Calculate S-G coefficients

HalfWin = ((F+1)/2) -1;

for n = (F+1)/2:(length(x)-5)-(F+1)/2,

% Zeroth derivative (smoothing only)
  SG0(n) = dot(g(:,1),y(n - HalfWin:n + HalfWin));

end

plot(x(1:length(SG0)),[y(1:length(SG0))',SG0']); % plot of noisy and smoothed signal

axis ([480 520 -20 20]); % changes view of axis to certain limits

set(gca,'Color',[0.8 0.8 0.8]); % changes background colour to grey

ylabel('Noisy Signal')

Questions:

- How can I change the colour and width of the lines?

(I have tried using the code below but this changes both line colours) >>plot(x(1:length(SG0)),[y(1:length(SG0))',SG0'],'k','LineWidth',2.5);

- How can I present other filters (Different N or F value on the same figure)

(I am able to produce multiple smoothed filters on top of the noisy signal but I'd want it to be already in colour within the codes rather than manually editing it in)

- Rather than visually seeing which smoothing filter best represents the signal. Is there another way to compare the two signals (noisy + smoothed) with one another and give a numerical value at the end?

(I tried to use corrcoef (x,y) but I do not know how to apply it properly)

I hope I can get some help from you guys.

Thanks

Answer 1

[Mohammad Abouali]

do it like this

plot(x(1:length(SG0)),y(1:length(SG0))','k','LineWidth',2.5);
hold on
plot(x(1:length(SG0)),SG0','r','LineWidth',2.5);

you can add more plots after hold on

plot(x(1:length(SG0)),FilteredDataUsingOtherParameters,'r','LineWidth',2.5);

then you can add the legend

legend('noisy data','Filtered using parameter 1', 'Filtered Using parameter 2')

corrcoef gives correlation coefficient. Once you smooth your data the correlation coefficient goes down anyway. The more it is smoothed the lower it would be. So, it is kinda hard to decide what is good and what is bad. Because the one that has the highest correlation coefficient is the one that is not smoothed at all (i.e. your original data), which certainly is the one that you don't want.

[Vince Roman]

Hi,

Thank you very much, I have managed to produce the plot of different filters on the same noisy signal. However I think the way I am running the codes I am over writing the previous smoothing filter and just plotting the results. Will this affect the comparison when I choose to compare the noisy signal with any of the filters?

This is how I currently run the codes:

load 'SineNoise'

Filter1

Filter2

Filter3

Filter4

Filter5

Plot

1) load 'SineNoise'

- This is the same from previous post

2) Filter1

N = 2; % Order of polynomial fit 2,4,6,8

F = 11; % Window length 11,23,35,47,59

[b,g] = sgolay(N,F); % Calculate S-G coefficients

HalfWin = ((F+1)/2) -1;

for n = (F+1)/2:(length(x)-5)-(F+1)/2,

% Zeroth derivative (smoothing only) SG0(n) = dot(g(:,1),y(n - HalfWin:n + HalfWin)); end

%Plot against x 0,0.1,0.2,....,

plot(x(1:length(SG0)),y(1:length(SG0)));

hold on

plot(x(1:length(SG0)),SG0','y','LineWidth',2.5);

3) Filter2 to Filter5 are pretty much the same except for the F value changing.

N = 2; % Order of polynomial fit 2,4,6,8

F = 23; % Window length 11,23,35,47,59

[b,g] = sgolay(N,F); % Calculate S-G coefficients

HalfWin = ((F+1)/2) -1;

for n = (F+1)/2:(length(x)-5)-(F+1)/2,

% Zeroth derivative (smoothing only) SG0(n) = dot(g(:,1),y(n - HalfWin:n + HalfWin)); end

%Plot against x 0,0.1,0.2,....,

plot(x(1:length(SG0)),SG0','m','LineWidth',2.5);

4) Plot - changing the axis and background colour aswell as legend

axis ([480 520 -20 20]); %changes view of axis to certain limits

legend('Noisy data','Filter1', 'Filter2', 'Filter3', 'Filter4', 'Filter5')

set(gca,'Color',[0.8 0.8 0.8]); %changes background colour to grey

ylabel('Noisy Signal')

I am just curious if it's possible when comparing the filter with the noisy signal to do something like this: Difference between the filter and the noise for each data point at a certain range (say from 480-520) and the average between the two. I am basically just looking for a number or a percentage difference between the two waveforms.

Or is there another/better way to compare the waveforms?

Thanks

[Image Analyst]

Overwriting a previous version of a signal does not affect the comparison of the current signal with your original noisy waveform - why would it?

Your proposed comparison, the "mean absolute difference" (see Wikipedia) is fine.