Savitzky-Golay algoritmus
A Savitzky-Golay algoritmus a jelfeldolgozás során alkalmazott módszer a görbe simítására és az egymást követő származékok kivonására . Abraham Savitzky (en) és Marcel Golay írták le 1964- ben .
Az algoritmus leírása
Tekintsünk egy y = ƒ ( x ) görbét , amely „asperitásokat”, alacsony amplitúdójú oszcillációkat mutat be; zajos jelzésről beszélünk . Ez egy diszkrét görbe, azaz pontok felhője ( x ( i ), y ( i )) 1 ≤ i ≤ n határozza meg .
A legegyszerűbb simítási algoritmus a futóátlagok módszere :
- ablakot, intervallumot tekintünk „félszélességűnek” ℓ (pontok száma);
- kiszámítjuk a függvény y ℓ + 1 átlagát az [1; 2 × ℓ + 1] (az intervallum szélessége tehát 2 ℓ + 1, ℓ nem pontosan a félszélesség);
- definiálunk egy új pontot ( x ( ℓ + 1), y liss ( ℓ + 1)) y liss ( ℓ + 1) = y ℓ + 1 ;
- az intervallumot egy ponttal húzzuk, és kezdjük újra.
Ez az algoritmus egyszerű, de „erőszakos”: nagymértékben csillapítja az erős variációkat, elcsípi a csúcsokat.
A következőkben az i intervallumot [ i - ℓ ; i + ℓ ], az i középső intervallum és a szélesség 2 ℓ + 1.
Kifinomultabb módon áll megfontolni egy polinom foka d , ahol d <2 ℓ + 1 minden intervallum i végzünk a regresszió , hogy meghatározza a polinom P i minimalizálja a hiba abban az értelemben, a legkisebb négyzetek. Meghatározzuk a simított értéket
ysima=Pén(x(én)).{\ displaystyle y _ {\ textrm {liss}} = P_ {i} (x (i)).}Sőt, ha a polinom legalább 1 fokos, akkor meghatározhatjuk a deriváltat
ysima′=Pén′(x(én)).{\ displaystyle y '_ {\ textrm {liss}} = P_ {i}' (x (i)).}és általában meghatározhatjuk a d- deriváltat.
Látjuk, hogy a csúszó eszköz módszer Savitzky-Golay simítás 0 fokú polinommal.
A paraméterek megválasztása
Gyakorlatban :
- a 2. fokú polinom lehetővé teszi a görbület figyelembevételét ; a 3. fokú polinom lehetővé teszi az inflexiós pontok figyelembevételét ; ritkán kell túllépni;
- az intervallumban lévő pontok számának elég nagynak kell lennie a polinom mértékéhez képest, hogy a simítás hatékony legyen; ha egyenlőek, akkor a polinom pontosan áthalad az intervallum minden pontján, tehát nincs simítás.
Ezért általában 3 fokú polinomot és legalább 5 pontos csúszó ablakot veszünk fel.
Minél nagyobb az ablak, annál simább a görbe, de annál kisebb a hullámhossz-variáció. A határszélesség megegyezik a megtartani kívánt részletek hullámhosszának sorrendjével. Ha például csúcsok érdekelnek minket - nagyon gyakori esetek a spektrometriában és a diffraktometriában -, akkor az ablak szélességét általában a legvékonyabb csúcs középmagasságának szélességének vesszük.
Állandó lépés esete
Konkrétan a polinom regresszió együtthatóit a csúszó intervallum pontjainak lineáris kombinációjával számoljuk . A simított és a derivált értékek tehát önmagukban lineáris kombinációk:
ylénss én=b0⋅yén-l+b1⋅yén-l+1+...+bl⋅yén+...+b2l⋅yén+l=∑k=02lbk⋅yén-l+k{\ displaystyle y _ {\ mathrm {liss} \ i} = b_ {0} \ cdot y_ {il} + b_ {1} \ cdot y_ {i-l + 1} + \ ldots + b_ {l} \ cdot y_ { i} + \ ldots + b_ {2l} \ cdot y_ {i + l} = \ sum _ {k = 0} ^ {2l} b_ {k} \ cdot y_ {il + k}} ;
ylénss én′=b0′⋅yén-l+b1′⋅yén-l+1+...+bl′⋅yén+...+b2l′⋅yén+l=∑k=02lbk′⋅yén-l+k{\ displaystyle y '_ {\ mathrm {liss} \ i} = b' _ {0} \ cdot y_ {il} + b '_ {1} \ cdot y_ {i-l + 1} + \ ldots + b' _ {l} \ cdot y_ {i} + \ ldots + b '_ {2l} \ cdot y_ {i + l} = \ sum _ {k = 0} ^ {2l} b' _ {k} \ cdot y_ { il + k}} ;
...
Ha a h = x i - x i - 1 lépés állandó, akkor a b k , b ' k együtthatók mindenhol azonosak; véges impulzusválaszűrő (FIR szűrő) esetén találjuk magunkat . Ezért kezdetben meghatározhatjuk a b k , b ' k ,… együtthatókat , amelyek gyors és könnyen megvalósítható algoritmust biztosítanak - például egy táblázattal ; ezután konvolúciós együtthatókról beszélünk .
Ilyen feltételek mellett azt tapasztaljuk, hogy:
- a simítás megegyezik, függetlenül attól, hogy 2-es (kvadratikus) vagy 3-as (köbös) fokú polinomot használunk-e;
- a levezetés ugyanaz, függetlenül attól, hogy 3 (köbös) vagy 4 („kvartikus”) fokú polinomot használunk;
- a második derivált ugyanaz, függetlenül attól, hogy 4 ("kvartikus") vagy 5 ("kvintikus") fokú polinomot használunk-e;
- az együtthatók szimmetrikusak vagy antiszimmetrikusak, a levezetés sorrendjének megfelelően:
b k = b 2 ℓ - k ; b ' k = - b' 2 ℓ - k ; ...
- mindig b ' l = 0: az x i pont nem járul hozzá y' liss i kiszámításához .
Az együtthatók értékei táblázatosak az adott ablakméretekhez (csúszó intervallum). Például egy 5 pontos ablak és egy 3 fokú polinom esetén:
ylénss én=135(-3⋅yén-2+12.⋅yén-1+17.⋅yén+12.⋅yén+1-3⋅yén+2){\ displaystyle y _ {\ mathrm {liss} \ i} = {\ frac {1} {35}} (- 3 \ cdot y_ {i-2} +12 \ cdot y_ {i-1} +17 \ cdot y_ {i} +12 \ cdot y_ {i + 1} -3 \ cdot y_ {i + 2})} ;
ylénss én′=112.h(1⋅yén-2-8.⋅yén-1+0⋅yén+8.⋅yén+1-1⋅yén+2){\ displaystyle y '_ {\ mathrm {liss} \ i} = {\ frac {1} {12h}} (1 \ cdot y_ {i-2} -8 \ cdot y_ {i-1} +0 \ cdot y_ {i} +8 \ cdot y_ {i + 1} -1 \ cdot y_ {i + 2})} ;
ylénss én′′=17h2(2⋅yén-2-1⋅yén-1-2⋅yén-1⋅yén+1+2⋅yén+2){\ displaystyle y _ {\ mathrm {liss} \ i} ^ {\ prime \ prime} = {\ frac {1} {7h ^ {2}}} (2 \ cdot y_ {i-2} -1 \ cdot y_ {i-1} -2 \ cdot y_ {i} -1 \ cdot y_ {i + 1} +2 \ cdot y_ {i + 2})}.
Az alábbi táblázatok megadják a konvolúciós együtthatókat adott polinomi fokokra és ablakszélességekre.
Konvolúciós együtthatók a simításhoz
Fokozat
|
2/3 (másodfokú / köbös)
|
4/5 (kvartikus / kvintikus)
|
---|
Szélesség
|
5. |
7 |
9.
|
7 |
9.
|
---|
–4
|
|
|
–21 |
|
15
|
---|
–3
|
|
–2 |
14 |
5. |
–55
|
---|
–2
|
–3 |
3 |
39 |
-30 |
30
|
---|
–1
|
12. |
6. |
54. |
75 |
135
|
---|
0
|
17. |
7 |
59 |
131 |
179
|
---|
1
|
12. |
6. |
54. |
75 |
135
|
---|
2
|
–3 |
3 |
39 |
-30 |
30
|
---|
3
|
|
–2 |
14 |
5. |
–55
|
---|
4
|
|
|
–21 |
|
15
|
---|
Szabványosítás
|
35 |
21 |
231 |
231 |
429
|
---|
A derivált konverziós együtthatói
Fokozat
|
2 (másodfokú)
|
3/4 (köbös / kvartikus)
|
---|
Szélesség
|
5. |
7 |
9.
|
5. |
7 |
9.
|
---|
–4
|
|
|
–4 |
|
|
86
|
---|
–3
|
|
–3 |
–3 |
|
22. |
–142
|
---|
–2
|
–2 |
–2 |
–2 |
1 |
–67 |
–193
|
---|
–1
|
–1 |
–1 |
–1 |
–8 |
–58 |
–126
|
---|
0
|
0 |
0 |
0 |
0 |
0 |
0
|
---|
1
|
1 |
1 |
1 |
8. |
58 |
126.
|
---|
2
|
2 |
2 |
2 |
–1 |
67 |
193
|
---|
3
|
|
3 |
3 |
|
–22 |
142
|
---|
4
|
|
|
4 |
|
|
–86
|
---|
Szabványosítás
|
10h |
28h |
60h |
12h |
252h |
1 188h
|
---|
Megjegyzés: a normalizálás a h mintavételi lépést használja.
Konvolúciós együtthatók a második deriváltra
Fokozat
|
2/3 (másodfokú / köbös)
|
4/5 (kvartikus / kvintikus)
|
---|
Szélesség
|
5. |
7 |
9.
|
5. |
7 |
9.
|
---|
–4
|
|
|
28. |
|
|
−4 158
|
---|
–3
|
|
5. |
7 |
|
–117 |
12,243
|
---|
–2
|
2 |
0 |
–8 |
–3 |
603 |
4,983
|
---|
–1
|
–1 |
–3 |
–17 |
48 |
–171 |
−6 963
|
---|
0
|
–2 |
–4 |
–20 |
–90 |
–630 |
−12 210
|
---|
1
|
–1 |
–3 |
–17 |
48 |
–171 |
−6 963
|
---|
2
|
2 |
0 |
–8 |
–3 |
603 |
4,983
|
---|
3
|
|
5. |
7 |
|
–117 |
12,243
|
---|
4
|
|
|
28. |
|
|
−4 158
|
---|
Szabványosítás
|
7 óra 2 |
42 óra 2 |
462 óra 2 |
36 óra 2 |
1 188 óra 2 |
56 628 óra 2 |
---|
A konvolúciós együtthatók kiszámítása
Egy adott k ponthoz megyünk . A kísérleti pont koordinátái ( x k , y k ).
A korrelációs együtthatók meghatározásához változó változást hajtunk végre annak érdekében, hogy csökkentett középpontú abszcisszája legyen - vagyis az az intervallum, amely alatt a regressziót elvégezzük, [- ℓ ; ℓ ]: pózolunk
z=x-x¯h{\ displaystyle z = {\ frac {x - {\ bar {x}}} {h}}}vagy
-
x az átlagos értékeinek x felett intervallum [ x k - ℓ ; x k + ℓ ] figyelembe véve, ez az ablak közepe;
-
h a mintavételi lépés (két szomszédos pont távolsága).
A z abszcissza az (- ℓ ; - ℓ + 1;…; –1; 0; 1;… ℓ ) értékeket veszi fel . Például egy 5 pontos csúszóablak esetén ( z i ) 1 ≤ i ≤ 5 = (–2; –1; 0; 1; 2).
A d fok polinomját ekkor fejezzük ki:
y=P(z)=nál nél0+nál nél1⋅z+nál nél2⋅z2+...+nál néld⋅zd{\ displaystyle y = P (z) = a_ {0} + a_ {1} \ cdot z + a_ {2} \ cdot z ^ {2} + \ ldots + a_ {d} \ cdot z ^ {d}}Tekintsük a vektor függvényt
F:(nál nél1nál nél2⋮nál néld)⟼(P(z1)P(z2)⋮P(z2l+1)){\ displaystyle \ mathrm {F}: {\ begin {pmatrix} a_ {1} \\ a_ {2} \\\ vdots \\ a_ {d} \ end {pmatrix}} \ longmapsto {\ begin {pmatrix} \ mathrm {P} (z_ {1}) \\\ mathrm {P} (z_ {2}) \\\ vdots \\ P (z_ {2l + 1}) \ end {pmatrix}}}val vel
(P(z1)P(z2)⋮P(z2l+1))≃(yk-lyk-l+1⋮yk+l){\ displaystyle {\ begin {pmatrix} P (z_ {1}) \\ P (z_ {2}) \\\ vdots \\ P (z_ {2l + 1}) \ end {pmatrix}} \ simeq {\ kezdet {pmatrix} y_ {kl} \\ y_ {kl + 1} \\\ vdots \\ y_ {k + l} \ end {pmatrix}}}A Jacobi mátrix írva
J=(∂P∂nál nélj(zén))=(1z1z12⋯z1d1z2z22⋯z2d⋮⋮⋮⋮1z2l+1z2l+12⋯z2l+1d){\ displaystyle \ mathrm {J} = \ balra ({\ frac {\ részleges P} {\ részleges a_ {j}}} (z_ {i}) \ jobbra) = {\ kezdődik {pmatrix} 1 & z_ {1 } & z_ {1} ^ {2} & \ cdots & z_ {1} ^ {d} \\ 1 & z_ {2} & z_ {2} ^ {2} & \ cdots & z_ {2} ^ {d } \\\ vdots & \ vdots & \ vdots && \ vdots \\ 1 & z_ {2l + 1} & z_ {2l + 1} ^ {2} & \ cdots & z_ {2l + 1} ^ {d} \ vége {pmatrix}}}A normál egyenletek megoldása
tJJ(nál nélén)=tJ(yj){\ displaystyle {} ^ {\ mathrm {t}} \ mathrm {J} \ mathrm {J} (a_ {i}) = {} ^ {\ mathrm {t}} \ mathrm {J} (y_ {j} )}ezért
(nál nélén)=(tJJ)-1tJ(yj){\ displaystyle (a_ {i}) = \ balra ({} ^ {\ mathrm {t}} \ mathrm {J} \ mathrm {J} \ jobbra) ^ {- 1} {} ^ {\ mathrm {t} } \ mathrm {J} (y_ {j})}A simított görbe értéke a polinom értéke az intervallum közepén:
ysima=P(0)=nál nél0{\ displaystyle y _ {\ textrm {liss}} = P (0) = a_ {0}} ;
Az egymást követő származékok esetében:
ylénss′=∂P∂x(xk)=∂z∂x(0)∂P∂z(0)=1hnál nél1{\ displaystyle y '_ {\ mathrm {liss}} = {\ frac {\ részleges P} {\ részleges x}} (x_ {k}) = {\ frac {\ részleges z} {\ részleges x}} ( 0) {\ frac {\ részleges P} {\ részleges z}} (0) = {\ frac {1} {h}} a_ {1}} ;
ylénss′′=∂2P∂x2(xk)=∂2z∂x2(0)∂P∂z(0)+∂z∂x∂2P∂z2(0)=2h2nál nél2{\ displaystyle y _ {\ mathrm {liss}} ^ {\ prime \ prime} = {\ frac {\ részleges ^ {2} \ mathrm {P}} {\ részleges x ^ {2}}} (x_ {k }) = {\ frac {\ részleges ^ {2} z} {\ részleges x ^ {2}}} (0) {\ frac {\ részleges \ mathrm {P}} {\ részleges z}} (0) + {\ frac {\ részben z} {\ részleges x}} {\ frac {\ részleges ^ {2} \ mathrm {P}} {\ részleges z ^ {2}}} (0) = {\ frac {2} {h ^ {2}}} a_ {2}}.
Például d = 3 fokú polinom és 5 pontos ablak esetén tehát therefore = 2:
J=(1-24-8.1-11-1100011111248.) ; tJ=(11111-2-101241014-8.-1018.){\ displaystyle \ mathrm {J} = {\ kezdődik {pmatrix} 1 & -2 & 4 & -8 \\ 1 & -1 & 1 & -1 \\ 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 1 & 2 & 4 & 8 \ end {pmatrix}} {\ text {; }} {} ^ {\ mathrm {t}} \ mathrm {J} = {\ kezdődik {pmatrix} 1 & 1 & 1 & 1 & 1 \\ - 2 & -1 & 0 & 1 & 2 \\ 4 & 1 & 0 & 1 & 4 \\ - 8 & -1 & 0 & 1 & 8 \ end {pmatrix}}}
(nál nél0nál nél1nál nél2nál nél3)=(tJJ)-1tJ(yk-2yk-1ykyk+1yk+2)=(-3/3512./3517./3512./35-3/351/12.-8./12.08./12.-1/12.2/14-1/14-2/14-1/142/14-1/12.2/12.0-2/12.1/12.)(yk-2yk-1ykyk+1yk+2){\ displaystyle {\ begin {pmatrix} a_ {0} \\ a_ {1} \\ a_ {2} \\ a_ {3} \ end {pmatrix}} = \ left ({} ^ {\ mathrm {t} } \ mathrm {J} \ mathrm {J} \ right) ^ {- 1} {} ^ {\ mathrm {t}} \ mathrm {J} {\ begin {pmatrix} y_ {k-2} \\ y_ { k-1} \\ y_ {k} \\ y_ {k + 1} \\ y_ {k + 2} \ end {pmatrix}} = {\ begin {pmatrix} -3 / 35 & 12/35 & 17 / 35 & 12/35 & - 3/35 \\ 1/12 & -8 / 12 & 0 & 8/12 & -1 / 12 \\ 2/14 & -1 / 14 & -2 / 14 & -1 / 14 & 2/14 \\ - 1/12 & 2/12 & 0 & -2 / 12 & 1/12 \ end {pmatrix}} {\ begin {pmatrix} y_ {k-2} \\ y_ {k- 1} \\ y_ {k} \\ y_ {k + 1} \\ y_ {k + 2} \ end {pmatrix}}}.
Tehát van:
ylénss k=nál nél0=135(-3yk-2+12.yk-1+17.yk+12.yk+1-3yk+2){\ displaystyle y _ {\ mathrm {liss} \ k} = a_ {0} = {\ frac {1} {35}} (- 3y_ {k-2} + 12y_ {k-1} + 17y_ {k} + 12év_ {k + 1} -3év_ {k + 2})} ;
ylénss k′=nál nél1h=112.h(yk-2-8.yk-1+0yk+8.k+1-yk+2){\ displaystyle y '_ {\ mathrm {liss} \ k} = {\ frac {a_ {1}} {h}} = {\ frac {1} {12h}} (y_ {k-2} -8y_ { k-1} + 0y_ {k} + 8_ {k + 1} -y_ {k + 2})} ;
ylénss k′′=2h2nál nél2=17h2(2yk-2-yk-1-2yk-yk+1+2yk+2){\ displaystyle y _ {\ mathrm {liss} \ k} ^ {\ prime \ prime} = {\ frac {2} {h ^ {2}}} a_ {2} = {\ frac {1} {7h ^ {2}}} (2y_ {k-2} -y_ {k-1} -2y_ {k} -y_ {k + 1} + 2y_ {k + 2})}.
Megjegyzések és hivatkozások
-
Abraham Savitzky és Marcel JE Golay , „ Az adatok simítása és differenciálása egyszerűsített legkisebb négyzetes eljárásokkal ”, Analitical Chemistry , vol. 8, n o 36,1964, P. 1627–1639 ( DOI 10.1021 / ac60214a047 )
-
nem mindig van állandó lépésünk, lehet változó lépésünk
- optimalizálási okokból: az egyik területet finomabban, a másikat durvábban akarjuk mérni (időmegtakarítás céljából);
- lényeges okokból:
- az x bemeneti értéket előállító rendszernek nincs lineáris válasza,
- nagyon precíz mérésekhez a szabályozási hurok hosszú időt vehet igénybe a bemeneti érték stabilizálásához az x th alapjelnél , ezután időt takaríthatunk meg egy olyan pont mérésével, amikor az alapérték nem érhető el teljesen, az x értékét nagyon pontosan rögzítjük.
<img src="https://fr.wikipedia.org/wiki/Special:CentralAutoLogin/start?type=1x1" alt="" title="" width="1" height="1" style="border: none; position: absolute;">