How can I find the filter bandwidth of Savitzky-Golay filtering?

The matlab function sgolayfilt(x,order,framelength) takes polynomial order and framelength as input arguments. As savitzky acts as a lowpass filter. Is there a possibility to find the filter bandwidth based on the input arguments (order and framelength)?

Answer 1

[Bruno Luong]

sa-golay filter is just non-causal FIR filter, where the coefficients are computed as following

order=3;framelength=11;
B = sgolay(order,framelength,[]);
b = B((framelength-1)./2+1,:);
a = 1;

Here is the frequency response using freqz function

[h,w] = freqz(b,a,2001);
plot(w/pi,20*log10(abs(h)))
ax = gca;
ax.YLim = [-100 20];
ax.XTick = 0:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

[Bjorn Gustavsson]

That was a good reminder for me! This is one of the good things with participating here, I regularly get reminded about functions and functionalities that has receded into the shadowy parts of memory. I completely forgot about that function (I admit to having an irrational and unreasonable aversion to SG-filtering, I should know better).

[Kalasagarreddi Kottakota]

Hi @Bruno Luong, Thanks for the answer. Can the filter bandwidth be found automatically using the bandwidth() function?

[Bruno Luong]

I don't know bandwidth function

[Star Strider]

The bandwidth function is a Control Systems and System Identification Toolbox function.

order=3;framelength=11;
B = sgolay(order,framelength,[]);
b = B((framelength-1)./2+1,:);
a = 1;

sgf = tf(b,a,1);                                        % Discrete Filter, Ts=1

BW = bandwidth(sgf)

BW = 0.6192

figure
bodeplot(sgf, {0 3})
grid

.

[Paul]

I realize that bandwidth is focused on the gain curve, but just want to point out that both solutions in this thread need to be adjusted to account for the fact that the SG filter is noncausal IF the phase response is also desired.

order=3;framelength=11;
B = sgolay(order,framelength,[]);
b = B((framelength-1)./2+1,:);
a = 1;
[h,w] = freqz(b,a,2001);
N = (framelength - 1)/2;
h = h.*exp(1j*w*N);
h = real(h);  % better

% best?
n = -N:N;
h = sum(b(:).*cos(n(:).*w.'));

sgf = tf(b,a,1)/tf('z',1)^N;  % or sgf = tf(b,a,-1)/tf('z',-1)^N

[Bjorn Gustavsson]

@Paul, is it even possible to get a causal filter using a window over a region centred such that

[Bruno Luong]

@Bjorn Gustavsson to me as soon as one use f(i+j) with j>0 it is non-causal.

This negative lag is translated by the phase shift of exp(1j*w*N) in the transfer function as Paul correctly expressed.

[Paul]

@Bjorn Gustavsson

I don't understand why you're asking about a causal filter. Did anything in my comment imply that the filter with coefficients b is causal?

To the contrary, the reason for the additional phase shift is that freqz assumes the inputs b,a describe a causal filter, but in this case it isn't. Hence the phase correction is needed.

Assume a simple case with a frame length of 5. The values in b are the coeffiecients of a noncausal filter of the form:

H(z) = b(1)*z^2 + b(2)*z + b(3) + b(4)*z^-1 + b(5)*z^-2

which can be rewritten as

H(z) = z^2 * (b(1) + b(2)*z^-1 + b(3)*z^-2 + b(4)*z^-3 + b(5)*z^-4)

With inputs b and a=1, freqz computes the response of the quantity inside the parentheses. The additional phase shift accounts for the z^2 outside the parentheses.

For completeness, the tf command, w/o additional arguments to specify otherwise, assumes that b, with a=1, describes a filter of the form

H(z) = b(1)*z^4 + b(2)*z^3 + b(3)*z^2 + b(4)*z + b(5)

which is why in that construction we'd divide by z^2.

Answer 2

[Bjorn Gustavsson]

QD suggestion:

x = randn(4096,1);
x_sgf = sgolayfilt(x,order,framelength);
fx = fftshift(fft(x));
fxsgf = fftshift(fft(x_sgf));
plot(medtilt1(abs(fxsgf)./abs(fx),3))

That should give you some sense of the filter characteristics. It is a rather brute-force, neither clever nor elegant way of calculating the frequency-response of the filter, but should give you a hint. The medfilt1 step I put in because the ratio became a bit spikey most likely due to numerical precision issues.

HTH

Answer 3

[Image Analyst]

I'm pretty sure the SG filter is a non-linear, non-causal filter like @Bruno Luong says. (The median filter is another such non-linear filter.) Thus there is no "frequency response" per se for the SG filter with certain parameters. You can only get the empirical frequency response if you know the actual input signal and compare actual, original signal to the filtered version of it. This is what @Bjorn Gustavsson showed. And of course it will vary depending on what your original input signal is.

[Bruno Luong]

I'm pretty sure the SG filter is a non-linear ...

Anything to backup your claim? I give the formula for the (linear) filter corresponds to SG. Ses B = sgolay(order,framelength,[]);

b = B((framelength-1)./2+1,:);

where B is pseudo-inverse of the Vandermonth matrix. b abd B are written as c and J in the wikipedia https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter (section derivation of convolustion coefficients). It is linear since it ths result of l^2 fitting of the polynomial of fix order through a fix number of data.

[Bjorn Gustavsson]

@Image Analyst, I'm of the understanding that the SG-filters are linear filters (at least I recall reading it that way and getting a proper "aha-moment" when I read that the polynomial fitting actually corresponded to a normal convolution with clever filtering kernels), are you sure that you're right about the non-linear part?