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