Фильтр Гаусса

* Тема: Реализация гауссова фильтра.
* Ключевые слова: фильтрация, фильтр, сглаживание, амплитудно-частотная характеристика, Гаусс, гауссовский, гауссовый, гауссов.
* Автор: Рейналь Левек.
* Перевод: А. Балабанов.
* Опубликован: 22.03.2008.
* Размещён: http://www.spsstools.ru/Syntax/TimeSeries/GaussianFilter.sps.txt (.sps).
* Проверен: SPSS 15.0.1.1.

* Пример гауссовского фильтра.
* По материалам: Notes_8, GEOS 595E, Spring 2003.
* (см. аналогичные материалы http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf ;
   из этого документа будут ясны некоторые особенности алгоритма и константы,
   в нём использованные - А.Б.)
* См. также http://www.statistics.com/resources/glossary/g/gaussfilt.php - А.Б.

* Синтаксис написан Рейналем Левек 21.05.2004..
* Мой сайт, посвящённый работе с SPSS: www.spsstools.net.

* Использованы данные из статистического сборника.
* Среднегодовые урожаи пшеницы в Англии (1259-1938 гг.). Составлено Bureau of Crop Estimates (бюро сельскохозяйственных прогнозов).

DATA LIST LIST /Year(A10) 	Bushel(F6).
BEGIN DATA
12/31/1259	17
12/31/1260	14
12/31/1261	13
12/31/1262	18
12/31/1263	12
12/31/1264	13
12/31/1265	10
12/31/1266	13
12/31/1267	13
12/31/1268	16
12/31/1269	15
12/31/1270	19
12/31/1271	21
12/31/1272	19
12/31/1273	16
12/31/1274	20
12/31/1275	15
12/31/1276	19
12/31/1277	15
12/31/1278	13
12/31/1279	15
12/31/1280	15
12/31/1281	18
12/31/1282	18
12/31/1283	21
12/31/1284	15
12/31/1285	16
12/31/1286	14
12/31/1287	9
12/31/1288	9
12/31/1289	13
12/31/1290	19
12/31/1291	17
12/31/1292	16
12/31/1293	25
12/31/1294	27
12/31/1295	20
12/31/1296	14
12/31/1297	16
12/31/1298	16
12/31/1299	18
12/31/1300	14
12/31/1301	15
12/31/1302	15
12/31/1303	12
12/31/1304	17
12/31/1305	15
12/31/1306	12
12/31/1307	17
12/31/1308	21
12/31/1309	23
12/31/1310	21
12/31/1311	13
12/31/1312	15
12/31/1313	17
12/31/1314	25
12/31/1315	45
12/31/1316	48
12/31/1317	25
12/31/1318	14
12/31/1319	17
12/31/1320	19
12/31/1321	35
12/31/1322	27
12/31/1323	22
12/31/1324	22
12/31/1325	17
12/31/1326	11
12/31/1327	12
12/31/1328	19
12/31/1329	20
12/31/1330	22
12/31/1331	24
12/31/1332	14
12/31/1333	13
12/31/1334	12
12/31/1335	16
12/31/1336	15
12/31/1337	11
12/31/1338	10
12/31/1339	18
12/31/1340	11
12/31/1341	11
12/31/1342	12
12/31/1343	17
12/31/1344	11
12/31/1345	11
12/31/1346	21
12/31/1347	20
12/31/1348	13
12/31/1349	16
12/31/1350	25
12/31/1351	31
12/31/1352	22
12/31/1353	13
12/31/1354	16
12/31/1355	18
12/31/1356	18
12/31/1357	21
12/31/1358	17
12/31/1359	18
12/31/1360	19
12/31/1361	16
12/31/1362	23
12/31/1363	26
12/31/1364	22
12/31/1365	18
12/31/1366	20
12/31/1367	26
12/31/1368	20
12/31/1369	36
12/31/1370	28
12/31/1371	21
12/31/1372	24
12/31/1373	19
12/31/1374	26
12/31/1375	23
12/31/1376	14
12/31/1377	11
12/31/1378	11
12/31/1379	17
12/31/1380	19
12/31/1381	17
12/31/1382	16
12/31/1383	15
12/31/1384	17
12/31/1385	15
12/31/1386	12
12/31/1387	10
12/31/1388	11
12/31/1389	16
12/31/1390	26
12/31/1391	16
12/31/1392	10
12/31/1393	11
12/31/1394	12
12/31/1395	15
12/31/1396	18
12/31/1397	17
12/31/1398	16
12/31/1399	17
12/31/1400	18
12/31/1401	22
12/31/1402	20
12/31/1403	15
12/31/1404	12
12/31/1405	11
12/31/1406	13
12/31/1407	14
12/31/1408	22
12/31/1409	27
12/31/1410	15
12/31/1411	15
12/31/1412	15
12/31/1413	13
12/31/1412	13
12/31/1415	19
12/31/1416	24
12/31/1417	16
12/31/1418	21
12/31/1419	15
12/31/1420	19
12/31/1421	16
12/31/1422	13
12/31/1423	13
12/31/1424	15
12/31/1425	12
12/31/1426	12
12/31/1427	13
12/31/1428	27
12/31/1429	24
12/31/1430	18
12/31/1431	14
12/31/1432	21
12/31/1433	18
12/31/1434	16
12/31/1435	17
12/31/1436	16
12/31/1437	28
12/31/1438	45
12/31/1439	23
12/31/1440	12
12/31/1441	12
12/31/1442	12
12/31/1443	13
12/31/1444	12
12/31/1445	19
12/31/1446	18
12/31/1447	16
12/31/1448	17
12/31/1449	16
12/31/1450	20
12/31/1451	20
12/31/1452	17
12/31/1453	15
12/31/1454	12
12/31/1455	16
12/31/1456	15
12/31/1457	17
12/31/1458	17
12/31/1459	15
12/31/1460	21
12/31/1461	22
12/31/1462	13
12/31/1463	11
12/31/1464	12
12/31/1465	14
12/31/1466	16
12/31/1467	16
12/31/1468	17
12/31/1469	19
12/31/1470	17
12/31/1471	17
12/31/1472	12
12/31/1473	12
12/31/1474	14
12/31/1475	16
12/31/1476	15
12/31/1477	20
12/31/1478	20
12/31/1479	18
12/31/1480	18
12/31/1481	26
12/31/1482	31
12/31/1483	22
12/31/1484	16
12/31/1485	14
12/31/1486	16
12/31/1487	17
12/31/1488	17
12/31/1489	18
12/31/1490	15
12/31/1491	20
12/31/1492	13
12/31/1493	12
12/31/1494	14
12/31/1495	12
12/31/1496	16
12/31/1497	15
12/31/1498	16
12/31/1499	14
12/31/1500	18
12/31/1501	25
12/31/1502	24
12/31/1503	19
12/31/1504	15
12/31/1505	15
12/31/1506	16
12/31/1507	17
12/31/1508	12
12/31/1509	9
12/31/1510	12
12/31/1511	17
12/31/1512	27
12/31/1513	18
12/31/1514	14
12/31/1515	20
12/31/1516	16
12/31/1517	19
12/31/1518	18
12/31/1519	22
12/31/1520	28
12/31/1521	23
12/31/1522	18
12/31/1523	17
12/31/1524	15
12/31/1525	16
12/31/1526	19
12/31/1527	41
12/31/1528	28
12/31/1529	29
12/31/1530	22
12/31/1531	25
12/31/1532	24
12/31/1533	23
12/31/1534	21
12/31/1535	31
12/31/1536	32
12/31/1537	21
12/31/1538	21
12/31/1539	17
12/31/1540	17
12/31/1541	27
12/31/1542	24
12/31/1543	28
12/31/1544	27
12/31/1545	47
12/31/1546	25
12/31/1547	15
12/31/1548	24
12/31/1549	49
12/31/1550	54
12/31/1551	71
12/31/1552	32
12/31/1553	30
12/31/1554	56
12/31/1555	66
12/31/1556	85
12/31/1557	25
12/31/1558	28
12/31/1559	33
12/31/1560	43
12/31/1561	47
12/31/1562	33
12/31/1563	59
12/31/1564	33
12/31/1565	32
12/31/1566	49
12/31/1567	33
12/31/1568	34
12/31/1569	35
12/31/1570	30
12/31/1571	37
12/31/1572	39
12/31/1573	79
12/31/1574	43
12/31/1575	48
12/31/1576	67
12/31/1577	61
12/31/1578	51
12/31/1579	39
12/31/1580	50
12/31/1581	49
12/31/1582	56
12/31/1583	52
12/31/1584	47
12/31/1585	63
12/31/1586	96
12/31/1587	84
12/31/1588	43
12/31/1589	59
12/31/1590	69
12/31/1591	61
12/31/1592	50
12/31/1593	55
12/31/1594	96
12/31/1595	116
12/31/1596	139
12/31/1597	171
12/31/1598	114
12/31/1599	71
12/31/1600	87
12/31/1601	80
12/31/1602	73
12/31/1603	80
12/31/1604	75
12/31/1605	78
12/31/1606	76
12/31/1607	87
12/31/1608	133
12/31/1609	121
12/31/1610	82
12/31/1611	91
12/31/1612	110
12/31/1613	118
12/31/1614	117
12/31/1615	102
12/31/1616	105
12/31/1617	119
12/31/1618	115
12/31/1619	90
12/31/1620	79
12/31/1621	77
12/31/1622	141
12/31/1623	137
12/31/1624	116
12/31/1625	124
12/31/1626	118
12/31/1627	92
12/31/1628	76
12/31/1629	104
12/31/1630	139
12/31/1631	157
12/31/1632	127
12/31/1633	132
12/31/1634	130
12/31/1635	126
12/31/1636	132
12/31/1637	130
12/31/1638	141
12/31/1639	107
12/31/1640	109
12/31/1641	121
12/31/1642	101
12/31/1643	121
12/31/1644	100
12/31/1645	116
12/31/1646	137
12/31/1647	159
12/31/1648	178
12/31/1649	188
12/31/1650	164
12/31/1651	148
12/31/1652	117
12/31/1653	78
12/31/1654	60
12/31/1655	84
12/31/1656	107
12/31/1657	117
12/31/1658	155
12/31/1659	159
12/31/1660	140
12/31/1661	165
12/31/1662	188
12/31/1663	130
12/31/1664	126
12/31/1665	115
12/31/1666	87
12/31/1667	89
12/31/1668	100
12/31/1669	116
12/31/1670	105
12/31/1671	102
12/31/1672	109
12/31/1673	132
12/31/1674	186
12/31/1675	135
12/31/1676	87
12/31/1677	106
12/31/1678	143
12/31/1679	134
12/31/1680	106
12/31/1681	117
12/31/1682	106
12/31/1683	106
12/31/1684	110
12/31/1685	110
12/31/1686	80
12/31/1687	86
12/31/1688	69
12/31/1689	85
12/31/1690	88
12/31/1691	90
12/31/1692	119
12/31/1693	169
12/31/1694	153
12/31/1695	132
12/31/1696	141
12/31/1697	162
12/31/1698	174
12/31/1699	169
12/31/1700	103
12/31/1701	85
12/31/1702	74
12/31/1703	85
12/31/1704	107
12/31/1705	81
12/31/1706	68
12/31/1707	69
12/31/1708	99
12/31/1709	199
12/31/1710	199
12/31/1711	137
12/31/1712	116
12/31/1713	122
12/31/1714	132
12/31/1715	108
12/31/1716	127
12/31/1717	112
12/31/1718	93
12/31/1719	93
12/31/1720	85
12/31/1721	89
12/31/1722	87
12/31/1723	85
12/31/1724	91
12/31/1725	113
12/31/1726	128
12/31/1727	110
12/31/1728	150
12/31/1729	138
12/31/1730	94
12/31/1731	72
12/31/1732	67
12/31/1733	70
12/31/1734	91
12/31/1735	107
12/31/1736	112
12/31/1737	107
12/31/1738	84
12/31/1739	95
12/31/1740	144
12/31/1741	130
12/31/1742	85
12/31/1743	63
12/31/1744	65
12/31/1745	66
12/31/1746	92
12/31/1747	92
12/31/1748	88
12/31/1749	60
12/31/1750	86
12/31/1751	99
12/31/1752	107
12/31/1753	113
12/31/1754	97
12/31/1755	80
12/31/1756	115
12/31/1757	181
12/31/1758	142
12/31/1759	104
12/31/1760	89
12/31/1761	73
12/31/1762	90
12/31/1763	99
12/31/1764	120
12/31/1765	139
12/31/1766	123
12/31/1767	190
12/31/1768	171
12/31/1769	133
12/31/1770	149
12/31/1771	152
12/31/1772	159
12/31/1773	160
12/31/1774	165
12/31/1775	152
12/31/1776	120
12/31/1777	143
12/31/1778	132
12/31/1779	105
12/31/1780	112
12/31/1781	140
12/31/1782	150
12/31/1783	165
12/31/1784	153
12/31/1785	131
12/31/1786	122
12/31/1787	129
12/31/1788	141
12/31/1789	160
12/31/1790	167
12/31/1791	148
12/31/1792	131
12/31/1793	150
12/31/1794	159
12/31/1795	229
12/31/1796	239
12/31/1797	163
12/31/1798	153
12/31/1799	210
12/31/1800	346
12/31/1801	363
12/31/1802	212
12/31/1803	179
12/31/1804	189
12/31/1805	273
12/31/1806	241
12/31/1807	229
12/31/1808	247
12/31/1809	296
12/31/1810	324
12/31/1811	290
12/31/1812	385
12/31/1813	334
12/31/1814	226
12/31/1815	199
12/31/1816	239
12/31/1817	295
12/31/1818	262
12/31/1819	227
12/31/1820	206
12/31/1821	171
12/31/1822	136
12/31/1823	162
12/31/1824	194
12/31/1825	208
12/31/1826	178
12/31/1827	178
12/31/1828	184
12/31/1829	202
12/31/1830	195
12/31/1831	202
12/31/1832	178
12/31/1833	161
12/31/1834	140
12/31/1835	120
12/31/1836	148
12/31/1837	170
12/31/1838	196
12/31/1839	215
12/31/1840	202
12/31/1841	196
12/31/1842	174
12/31/1843	152
12/31/1844	156
12/31/1845	155
12/31/1846	166
12/31/1847	212
12/31/1848	154
12/31/1849	135
12/31/1850	122
12/31/1851	117
12/31/1852	124
12/31/1853	162
12/31/1854	220
12/31/1855	227
12/31/1856	210
12/31/1857	171
12/31/1858	134
12/31/1859	133
12/31/1860	162
12/31/1861	168
12/31/1862	169
12/31/1863	136
12/31/1864	122
12/31/1865	127
12/31/1866	152
12/31/1867	196
12/31/1868	194
12/31/1869	146
12/31/1870	143
12/31/1871	172
12/31/1872	173
12/31/1873	178
12/31/1874	170
12/31/1875	137
12/31/1876	140
12/31/1877	173
12/31/1878	141
12/31/1879	133
12/31/1880	135
12/31/1881	138
12/31/1882	137
12/31/1883	126
12/31/1884	108
12/31/1885	100
12/31/1886	94
12/31/1887	99
12/31/1888	97
12/31/1889	90
12/31/1890	97
12/31/1891	113
12/31/1892	92
12/31/1893	80
12/31/1894	69
12/31/1895	70
12/31/1896	80
12/31/1897	92
12/31/1898	103
12/31/1899	78
12/31/1900	82
12/31/1901	81
12/31/1902	85
12/31/1903	81
12/31/1904	86
12/31/1905	90
12/31/1906	86
12/31/1907	93
12/31/1908	97
12/31/1909	111
12/31/1910	95
12/31/1911	95
12/31/1912	104
12/31/1913	95
12/31/1914	105
12/31/1915	159
12/31/1916	174
12/31/1917	230
12/31/1918	218
12/31/1919	219
12/31/1920	242
12/31/1921	216
12/31/1922	244
12/31/1923	127
12/31/1924	148
12/31/1925	157
12/31/1926	160
12/31/1927	148
12/31/1928	130
12/31/1929	127
12/31/1930	103
12/31/1931	74
12/31/1932	76
12/31/1933	69
12/31/1934	62
12/31/1935	67
12/31/1936	93
12/31/1937	120
12/31/1938	88
END DATA.
COMPUTE nobreak=1.
VAR LAB Year 'Год' / Bushel 'бушелей (млн?)'
SAVE OUTFILE='c:\\temp\\bushels.sav'.


************************************************************.
*** Задание параметров фильтра и вычисление вспомогательных переменных.
** Lambda05 - желаемая длина волны, для которой амплитудно-частотная характеристика фильтра составит 0.5.

* Возможно, уместно такое упрощенное понимание этого параметра: это тот отрезок времени,
  амплитуды колебаний внутри которого с длинами волн,
  равными этому отрезку времени, сократятся наполовину. А амплитуды ещё более коротких волн, разовых
  всплесков, сократятся ещё сильнее - А.Б.).

DEFINE !lambda05() 10 !ENDDEFINE.
DEFINE !var() bushel !ENDDEFINE.
************************************************************.

NEW FILE.
INPUT PROGRAM.
- COMPUTE #SigmaG=!lambda05 / 6.
- LOOP cnt=1 TO 101.
-   COMPUTE wgt=PDF.T((cnt-51)/#SigmaG,4300).
* 4300 - "достаточно большое" число степеней свободы, когда
  распределение Стьюдента практически совпадает со стандартным нормальным
  распределением. Аналогично можно было бы подсчитать веса так (примеч. пер.):
*   COMPUTE wgt=PDF.Norm((cnt-51)/#SigmaG,0,1).
-   END CASE.
- END LOOP.
- END FILE.
END INPUT PROGRAM.
FORMATS wgt (F9.4).
EXECUTE.
SAVE OUTFILE='c:\\temp\\weights.sav'.

* Найдём максимальное значение веса.

COMPUTE nobreak=1.
AGGREGATE OUTIFLE=*
	/PRESORTED
	/BREAK=nobreak
	/maxval=MAX(wgt).
WRITE OUTFILE='c:\\define max.sps' 
	/'DEFINE !maxval()' maxval(F8.7)'!ENDDEFINE.'.
EXECUTE.
INCLUDE 'c:\\define max.sps'.
GET FILE='c:\\temp\\weights.sav'.
* Отсеем "слишком малые" веса (примеч. пер.).
SELECT IF wgt>.05 * !maxval.
SUMMARIZE
  /TABLES=wgt
  /FORMAT=NOCASENUM VALID TOTAL
  /TITLE='Ненормированные веса'
  /MISSING=VARIABLE
  /CELLS=SUM .


* Вычислим сумму и количество оставшихся "крупных" весов.

COMPUTE nobreak=1.
AGGREGATE OUTIFLE=*
	/PRESORTED
	/BREAK=nobreak
	/sumval=SUM(wgt) /nbval=N.
COMPUTE nb=TRUNC(nbval/2).
WRITE OUTFILE='c:\\define sum.sps' 
	/'DEFINE !sumval()' sumval(F8.7)'!ENDDEFINE.'
	/'DEFINE !nbval()' nbval(F8)'!ENDDEFINE.'
	/'DEFINE !nb()' nb(F8)'!ENDDEFINE.'.
EXECUTE.
INCLUDE 'c:\\define sum.sps'.
title !sumval .
* Нормируем веса так, чтобы их сумма была равна 1.

GET FILE='c:\\temp\\weights.sav' /KEEP=wgt.
SELECT IF wgt>.05*!maxval.
COMPUTE wgt=wgt / !sumval.
*Проверим, что сумма равна 1.
SUMMARIZE
  /TABLES=wgt
  /FORMAT=NOCASENUM VALID TOTAL
  /TITLE='Сумма нормированных весов'
  /MISSING=VARIABLE
  /CELLS=SUM .
STRING w(a8).
COMPUTE w=CONCAT("w",LTRIM(STRING($CASENUM,F4))).
FLIP /NEWNAMES=w.

COMPUTE nobreak=1.
MATCH FILES FILE='c:\\temp\\bushels.sav'
	/TABLE=*
	/BY=nobreak
	/DROP=nobreak case_lbl.
SAVE OUTFILE='c:\\temp\\temp.sav'.

GET FILE='c:\\temp\\temp.sav'.



*///////////////.
DEFINE !doit (var=!TOKENS(1)  /nb=!TOKENS(1) /nbval=!TOKENS(1))

!LET !nbp2=!LENGTH(!CONCAT(!BLANKS(!nb),!BLANKS(2)))
!DO !cnt=1 !TO !nbval
  !IF (!BLANKS(!cnt) !LE !BLANKS(!nb)) !THEN
- COMPUTE !CONCAT('v',!cnt)=LAG(!var, !LENGTH(!SUBSTR(!BLANKS(!nb),!cnt)) ).
 !ELSE !IF ( !BLANKS(!cnt) !EQ !CONCAT(!BLANKS(!nb),!BLANKS(1)) ) !THEN
- COMPUTE !CONCAT('v',!cnt)=!var.
  !ELSE
- CREATE  !CONCAT('v',!cnt)=LEAD(!var, !LENGTH( !SUBSTR(!BLANKS(!cnt) , !nbp2)) ).
  !IFEND
  !IFEND
!DOEND
VECTOR v=v1 TO !CONCAT('v',!nbval) /w=w1 TO !CONCAT('w',!nbval).
NUMERIC gs.
VARIABLE LABEL gs !QUOTE(!CONCAT('Сглаживание фильтром Гаусса с параметром лямбда 0.5 =', !nbval))
!ENDDEFINE.
*///////////////.

EXECUTE.
SET MPRINT=yes.
!doit var=!var nb=!nb nbval=!nbval.

* Непосредственно вычисление сглаженных значений.

DO IF NVALID(v1 TO v9)=9.
- COMPUTE #gs=0.
- LOOP cnt=1 TO !nbval.
-   COMPUTE #gs=#gs + v(cnt)*w(cnt).
- END LOOP.
- COMPUTE gs=#gs.
END IF.

* Убираем лишние переменные.
ADD FILES FILE=* /KEEP=year bushel gs.
*Строим графики временных рядов. Сначала наложенные, затем - исходный, затем - сглаженный.

GRAPH
  /LINE(MULTIPLE)=MEAN(Bushel) MEAN(gs) BY Year
  /MISSING=LISTWISE .

TSPLOT VARIABLES= Bushel
  /NOLOG
  /FORMAT NOFILL NOREFERENCE.

TSPLOT VARIABLES= gs
  /NOLOG
  /FORMAT NOFILL NOREFERENCE.