1. 程式人生 > >Matlab中Savitzky-Golay filtering(最小二乘平滑濾波)函式sgolayfilt的使用方法

Matlab中Savitzky-Golay filtering(最小二乘平滑濾波)函式sgolayfilt的使用方法

語法規則

y = sgolayfilt(x,order,framelen)
y = sgolayfilt(x,order,framelen,weights)
y = sgolayfilt(x,order,framelen,weights,dim)

語法描述

y = sgolayfilt(x,order,framelen):對資料向量x使用Savitzky-Golay FIR平滑濾波器。如果x是一個矩陣,則sgolayfilt對每一列進行操作。多項式階次order必須小於框長度framelen,因此framelen必須是奇數。如果order=framelen-1,則濾波器不進行平滑。

y = sgolayfilt(x,order,framelen,weights):指定了一個權重向量weights,它為正實數,用於最小二乘估計的最小化。如果該值沒有指定,或定義為空[],它預設為單位矩陣。

y = sgolayfilt(x,order,framelen,weights,dim):指定了維數dim。如果dim沒有被指定,則sgolayfilt對第一個不為1的維進行操作;指定為1則對列向量進行操作,指定2則對行向量進行操作。

舉例

穩態和瞬態Savitzky-Golay濾波器

生成一個隨機訊號並使用sgolayfilt平滑。指定多項式階數為3,框長度為11.畫出原始的和平滑後的訊號。

order = 3;
framelen = 11;

lx = 34;
x = randn(lx,1);

sgf = sgolayfilt(x,order,framelen);

plot(x,':')
hold on
plot(sgf,'.-')
legend('signal','sgolay')

sgolayfilt函式通過對sgolay的輸出B的中間行進行卷積來執行大部分濾波。結果是過濾的訊號的穩態部分。生成並畫出這個部分。

m = (framelen-1)/2;

B = sgolay(order,framelen);

steady = conv(x,B(m+1,:),'same');

plot(steady)
legend('signal','sgolay','steady')

靠近訊號邊緣的樣本不能放置在對稱視窗的中心,並且必須以不同的方式進行處理。

為了確定啟動瞬態,矩陣將B的第一個(framelen-1)/2行乘以訊號的第一個framelen長度的取樣。

ybeg = B(1:m,:)*x(1:framelen);

為了確定終端瞬態,矩陣將B的最後的(framelen-1)/2行乘以訊號的最後的framelen長度的取樣。

yend = B(framelen-m+1:framelen,:)*x(lx-framelen+1:lx);

連線瞬態和穩態部分以產生完整的訊號。

cmplt = steady;
cmplt(1:m) = ybeg;
cmplt(lx-m+1:lx) = yend;

plot(cmplt)
legend('signal','sgolay','steady','complete')
hold off

新增權重會破壞B的對稱性,並且需要額外步驟來提供合適的解決方法。

對聲音訊號的Savitzky-Golay濾波

用三階Savitzky-Golay平滑mtlb訊號,資料框長度為41.

load mtlb
smtlb = sgolayfilt(mtlb,3,41);

subplot(2,1,1)
plot(1:2000, mtlb(1:2000))
axis([0 2000 -4 4])
title('mtlb')
grid

subplot(2,1,2)
plot(1:2000,smtlb(1:2000))
axis([0 2000 -4 4])
title('smtlb')
grid