* Тема: Реализация гауссова фильтра. * Ключевые слова: фильтрация, фильтр, сглаживание, амплитудно-частотная характеристика, Гаусс, гауссовский, гауссовый, гауссов. * Автор: Рейналь Левек. * Перевод: А. Балабанов. * Опубликован: 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.