sgolay

Diseño de filtros de Savitzky-Golay

Sintaxis

b = sgolay(m,fl)

b = sgolay(m,fl,w)

[b,g] = sgolay(___)

Descripción

b = sgolay(m,fl) diseña un filtro de suavizado FIR de Savitzky-Golay con orden polinomial m y longitud de cuadro fl.

b = sgolay(m,fl,w) especifica un vector de ponderación, w, que contiene las ponderaciones reales de valor positivo que se usarán durante la minimización de mínimos cuadrados.

[b,g] = sgolay(___) devuelve la matriz g de los filtros de diferenciación. Puede utilizar estos argumentos de salida con cualquiera de las sintaxis de entrada anteriores.

Ejemplos

Suavizado de Savitzky-Golay para sinusoide con ruido

Genere una señal que consiste en una sinusoide de 0,2 Hz integrada en ruido blanco gaussiano y muestreada cinco veces por segundo durante 200 segundos.

dt = 1/5;
t = (0:dt:200-dt)';

x = 5*sin(2*pi*0.2*t) + randn(size(t));

Utilice sgolay para suavizar la señal. Utilice marcos de 21 muestras y polinomios de cuarto orden.

order = 4;
framelen = 21;

b = sgolay(order,framelen);

Calcule la parte de estado estable de la señal convolucionándola con la fila central de b.

ycenter = conv(x,b((framelen+1)/2,:),'valid');

Calcule los transitorios. Utilice las últimas filas de b para el inicio y las primeras filas de b para el fin.

ybegin = b(end:-1:(framelen+3)/2,:) * x(framelen:-1:1);
yend = b((framelen-1)/2:-1:1,:) * x(end:-1:end-(framelen-1));

Concatene los transitorios y la parte de estado estable para generar la señal suavizada completa. Represente la señal original y la estimación de Savitzky-Golay.

y = [ybegin; ycenter; yend];
plot([x y])
legend('Noisy Sinusoid','S-G smoothed sinusoid')

Filtro de Savitzky-Golay con ponderaciones

Diseñe un filtro de suavizado de Savitzky-Golay con polinomios de quinto orden para filtrar una señal con ruido muestreada a ${\mathit{F}}_{\mathit{s}}=8192$ Hz. Utilice una ventana de Hamming. La longitud del cuadro es de 101 muestras.

Fs = 8192;
m = 5;
fl = 101;

weights = hamming(fl);

b = sgolay(m,fl,weights);

Represente la respuesta en frecuencia para la fila central de b. Utilice 1024 puntos de FFT.

NFFT = 1024;
freqz(b((fl+1)/2,:),1,NFFT,Fs)

Para obtener más información sobre los vectores de ponderación y su impacto en el rendimiento del filtrado de señales, consulte Design and Analyze Savitzky-Golay Filters.

Diferenciación de Savitzky-Golay

Genere una señal que consiste en una sinusoide de 0,2 Hz integrada en ruido blanco gaussiano y muestreada cuatro veces por segundo durante 20 segundos.

dt = 0.25;
t = (0:dt:20-1)';

x = 5*sin(2*pi*0.2*t)+0.5*randn(size(t));

Estime las tres primeras derivadas de la sinusoide utilizando el método de Savitzky-Golay. Utilice marcos de 25 muestras y polinomios de quinto orden. Divida las columnas por potencias de dt para escalar las derivadas correctamente.

[b,g] = sgolay(5,25);

dx = zeros(length(x),4);
for p = 0:3
  dx(:,p+1) = conv(x, factorial(p)/(-dt)^p * g(:,p+1), 'same');
end

Represente la señal original, la secuencia suavizada y las estimaciones de las derivadas.

plot(x,'.-')
hold on
plot(dx)
hold off

legend('x','x (smoothed)','x''','x''''', 'x''''''')
title('Savitzky-Golay Derivative Estimates')

Argumentos de entrada

m — Orden polinomial

Entero no negativo

Orden polinomial, especificado como entero no negativo. El valor de m debe ser menor que fl. Si m = fl – 1, el filtro diseñado no suaviza.

fl — Longitud del cuadro

entero positivo impar

Longitud del cuadro, especificada como entero impar positivo. El valor de fl debe ser mayor que m.

w — Vector de ponderación

ones(fl,1) (predeterminado) | vector real positivo

Vector de ponderación, especificado como un vector real positivo. El vector de ponderación tiene la misma longitud que fl y se utiliza para llevar a cabo la minimización de mínimos cuadrados.

Argumentos de salida

b — Coeficientes de filtro FIR variables en el tiempo

matriz

Coeficientes de filtro FIR variables en el tiempo, especificados como una matriz de fl por fl. En una implementación de filtro de suavizado (por ejemplo, sgolayfilt), las últimas filas (fl-1)/2 (cada una un filtro FIR) se aplican a la señal durante el transitorio de inicio, y las primeras filas (fl-1)/2 se aplican a la señal durante el transitorio de fin. La fila central se aplica a la señal en el estado estable.

g — Matriz de filtros de diferenciación

matriz

Matriz de filtros de diferenciación, especificada como una matriz. Cada columna de g es un filtro de diferenciación para derivadas de orden p-1, donde p es el índice de columnas. Dada una señal x de longitud fl, puede encontrar una estimación de la derivada de p.^-ésimo orden, xp, de su valor central de xp((fl+1)/2) = (factorial(p)) * g(:,p+1)' * x.

Algoritmos

Los filtros de Savitzky-Golay se utilizan para suavizar señales con ruido cuyo intervalo de frecuencia es grande. Los filtros de suavizado de Savitzky-Golay tienden a filtrar menos contenido de alta frecuencia de la señal que los filtros FIR promediadores estándar. Sin embargo, son menos eficaces para rechazar el ruido cuando los niveles de ruido son particularmente elevados.

En general, el filtrado consiste en sustituir cada punto de una señal por una combinación de los valores de señal contenidos en una ventana variable centrada en el punto, asumiendo que los puntos circundantes miden casi el mismo valor subyacente. Por ejemplo, los filtros de media móvil reemplazan cada punto de datos por la media local de los puntos de datos circundantes. Si un punto de datos dado tiene k puntos a la izquierda y k puntos a la derecha, para una longitud total de ventana de L = 2k + 1, el filtro de media móvil hace la sustitución

$$x_s\to {\widehat{x}}_s=\frac{1}{L}{\displaystyle \sum _{r=-k}^k{x}_{s+r}.}$$

Los filtros de Savitzky-Golay generalizan esta idea por mínimos cuadrados ajustando un polinomio de n-ésimo orden a través de los valores de señal en la ventana y tomando el punto central calculado de la curva del polinomio ajustado como el nuevo punto de datos suavizado. En un punto dado, x_s,

$$\left[\begin{array}{c}{x}_{s-k}\\ ⋮\\ x_{s-1}\\ x_s\\ x_{s+1}\\ ⋮\\ x_{s+k}\end{array}\right]=\left[\begin{array}{c}{b}_0+b_1\left(t_s-k\Delta t\right)+b_2{\left(t_s-k\Delta t\right)}^2+\cdots +b_n{\left(t_s-k\Delta t\right)}^n\\ ⋮\\ b_0+b_1\left(t_s-1\Delta t\right)+b_2{\left(t_s-1\Delta t\right)}^2+\cdots +b_n{\left(t_s-1\Delta t\right)}^n\\ b_0+b_1\left(t_s-0\Delta t\right)+b_2{\left(t_s-0\Delta t\right)}^2+\cdots +b_n{\left(t_s-0\Delta t\right)}^n\\ b_0+b_1\left(t_s+1\Delta t\right)+b_2{\left(t_s+1\Delta t\right)}^2+\cdots +b_n{\left(t_s+1\Delta t\right)}^n\\ ⋮\\ b_0+b_1\left(t_s+k\Delta t\right)+b_2{\left(t_s+k\Delta t\right)}^2+\cdots +b_n{\left(t_s+k\Delta t\right)}^n\end{array}\right]=\left[\begin{array}{c}{a}_0+a_1\left(-k\right)+a_2{\left(-k\right)}^2+\cdots +a_n{\left(-k\right)}^n\\ ⋮\\ a_0+a_1\left(-1\right)+a_2{\left(-1\right)}^2+\cdots +a_n{\left(-1\right)}^n\\ a_0+a_1\left(0\right)+a_2{\left(0\right)}^2+\cdots +a_n{\left(0\right)}^n\\ a_0+a_1\left(1\right)+a_2{\left(1\right)}^2+\cdots +a_n{\left(1\right)}^n\\ ⋮\\ a_0+a_1\left(k\right)+a_2{\left(k\right)}^2+\cdots +a_n{\left(k\right)}^n\end{array}\right]$$

o, en términos de matrices,

$$x=\left[\begin{array}{ccccc}1& -k& {\left(-k\right)}^2& \cdots & {\left(-k\right)}^n\\ ⋮& ⋮& ⋮& ⋰& ⋮\\ 1& -2& {\left(-2\right)}^2& \cdots & {\left(-2\right)}^n\\ 1& -1& {\left(-1\right)}^2& \cdots & {\left(-1\right)}^n\\ 1& 0& 0& \cdots & 0\\ 1& 1& 1^2& \cdots & 1^n\\ 1& 2& 2^2& \cdots & 2^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& k& k^2& \cdots & k^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ ⋮\\ a_n\end{array}\right]\equiv Ha.$$

Para encontrar las estimaciones de Savitzky-Golay, utilice la pseudoinversa de H para calcular a y luego multiplique previamente por H:

$$\widehat{x}=H{\left(H^{T}H\right)}^{-1}{H}^{T}x=Bx.$$

Para evitar condiciones inadecuadas, sgolay utiliza la función qr para calcular una descomposición de menor tamaño de H como H = QR, en cuanto a B = QQ^T.

Es necesario calcular B solo una vez. Las estimaciones de Savitzky-Golay en la mayoría de puntos de la señal son el resultado de convolucionar la señal con la fila central de B. El resultado es la parte de estado estable de la señal filtrada. Las primeras k filas de B producen el transitorio inicial y las últimas k filas de B producen el transitorio final. Para ver un ejemplo, consulte sgolayfilt. Es posible mejorar la supresión de ruido aumentando la longitud de la ventana, pero esto introduce un zumbido análogo al fenómeno de Gibbs cerca de los transitorios.

Referencias

[1] Orfanidis, Sophocles J. Introduction to Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1996.

[2] Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. New York: Cambridge University Press, 1992.

[3] Schafer, Ronald W. “What Is a Savitzky-Golay Filter? [Lecture Notes].” IEEE Signal Processing Magazine Vol. 28, Number 4, July 2011, pp. 111–117. https://doi.org/10.1109/MSP.2011.941097.

Capacidades ampliadas

Generación de código C/C++
Genere código C y C++ mediante MATLAB® Coder™.

Historial de versiones

Introducido antes de R2006a

Consulte también

Funciones

• fir1 | firls | filter | sgolayfilt

Temas

• Design and Analyze Savitzky-Golay Filters