авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 2 | 3 ||

«САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ На правах рукописи Усевич Константин Дмитриевич АНАЛИЗ ...»

-- [ Страница 4 ] --

Приведем расширение на двумерный случай способа ганкелизаци, пред ложенного в [81] для обычной ганкелизации. Определим (Nx + Lx 1) (Ny + Ly 1) массив следующим образом 0L 1Ly 1 0Lx 1Ky 0Lx 1Ly x Gj = 0Kx Ly 1.

j 0Kx,Ly 0Lx 1Ly 1 0Lx 1Ky 0Lx 1Ly 1.

Тогда Fj = (D) (Gj j ), где обозначает покомпонентное умножение, а элементы (D)mn = 1/#Dmn массива D MNx,Ny определены соотношениями (1.7.9) Ускорение построения ковариационной матрицы Рассмотрим сложность основных шагов алгоритма с построением ковари ационной матрицы (см. таблицу 5.8) Шаг количество умно жений O(L2 K) Вычисление ковариационной матрицы O(L3 ) Решение задачи собственных значений Вычисление p факторных векторов O(pLK) Вычисление r восстановленных массивов O(rLK) Таблица 5.8. Вычислительная сложность шагов алгоритма АСС В типичных задачах обработки изображений K = Kx Ky довольно вели ко. При выборе небольших размеров окна Lx Kx, Ly Ky, и небольшого числа компонент p, r Lx Ly для восстановления наиболее трудоемким ока зывается шаг построения ковариационной матрицы. Покажем, каким образом этот шаг можно ускорить, используя структуру ганкелевой матрицы. Отме тим, что для 0 u, s Lx 1 и 0 v, t Ly выполняется равенство s(u+1)+vLx +1,(s+1)+tLx +1 su+vLx +1,s+tLx +1 = (1,Ky ) (1,K ) (1,K ) = Fu,v y ), Fs,t (1,K + Fu+Kyx,v, Fs+Ky,t M, M x а для 0 u, s Lx и 0 v, t Ly 1 выполняется su+(v+1)Lx +1,s+(t+1)Lx +1 su+vLx +1,s+tLx +1 = (K,1) (K,1) (K,1) = Fu,vx,1), Fs,t x (K x x + Fu,v+Ky, Fs,t+Ky M.

M Таким образом, достаточно вычислить лишь всевозможные элементы s0,s+tLx +1, а остальные элементы smn можно вычислить итеративно, с помо щью 2 max(Kx, Ky ) умножений, что дает вычислительную сложность O(LK + L2 (Kx +Ky )) меньшую, чем O(L2 K). Отметим, что при данном способе вычис ления может накапливаться ошибка. Если Kx Ky, L Kx, то сложность вычисления ковариационной матрицы имеет порядок O(LK).

Заметим также, что ковариационная матрица может быть вычислена с по мощью процедуры быстрого вычисления произведения траекторной матрицы на вектор, описываемой в следующем разделе.

5.3.2. Быстрые вычисления с помощью БПФ В методах решения частичной задачи на собственные значения, основан ных на итерациях Ланцоша, построение траекторной матрицы не обязательно, и требуется лишь реализация операции умножения матрицы на вектор. В ра боте [81] была предложена схема вычисления АСС-разложения для рядов, в которой умножения матрицы на вектор эффективно вычисляются с помощью дискретного преобразования Фурье.

Дискретное преобразование Фурье и быстрое умножение ганкелевой матрицы на вектор Пусть N 1 и = e2i/N. Определим матрицу 1 1... 1 1... N FN =..

...

........

1 N 1... ( N 1 )N Дискретным преобразованием Фурье (ДПФ) комплексного вектора GN = def (g0,..., gN 1 )T называется вектор GN = FN GN. Матрица FN называется мат рицей ДПФ. Известно, что F1 = (FN ).

N N Для вектора GN обратным дискретным преобразованием Фурье будем назы вать вектор F1 GN. ДПФ и обратное ДПФ могут быть вычислены за O(N log N ) N умножений с помощью алгоритмов быстрого преобразования Фурье [78].

Пусть FN = (f0,..., fN 1 )T временной ряд и задана длина окна 1 L N, K = N L + 1. Определим вектор GN = (fK1,..., fN 1, f0,..., fK2 )T.

Тогда для матрицы GN = F1 diag(FN GN )FN (являющейся циркулянтной N матрицей [11] вектора GN ) и любых векторов U CL, V CK выполняются равенства X(L) (FN )V V =, GN 0L 0K =, GN (K) U (FN )U X где U и V перевернутые векторы U и V. Таким образом, для вычисле T ния произведения ганкелевых матриц X(L) (FN ) и X(K) (FN ) = X(L) (FN ) на вектор достаточно уметь вычислять произведение матрицы GN на вектор, которое можно выразить как GN A = F1 (FN A FN GN ), N что требует O(N log N ) операций. Подробнее см. в [81]. Заметим, что описан ный способ умножения является более экономичным по сравнению с умноже нием с помощью свертки, используемым в [104].

Двумерное преобразование Фурье Пусть G = G(Nx,Ny ) массив размера Nx Ny. Двумерное дискретное def преобразование массива (5.2.2) можно выразить как G = FNx GFT y, т.е. при N менение одномерного преобразования Фурье к строкам, а затем к столбцам.

Таким образом, двумерное дискретное преобразование Фурье можно вычис лить за Nx Ny log(Nx Ny ) операций.

По предложению 1.1.3, vec FNx GFNy = (FNy FNx ) vec G. Таким обра def зом, матрица F(Nx,Ny ) = FNy FNx является матрицей линейного оператора · ·.

Теперь рассмотрим массив F(Nx,Ny ) и размеры окна (Lx, Ly ), 1 Lx Nx, 1 Ly Ny, Kx = Nx Lx + 1, Ky = Ny Ly + 1.

Предложение 5.3.1. Пусть fKx 1,Ky 1... fKx 1,Ny 1 fKx 1,0... fKx 1,Ky....

........ def fNx 1,Ky 1... fNx 1,Ny 1... fNx 1,Ky fNx 1, G=.

f0,Ky 1... f0,Ny 1 f0,0... f0,Ky 2....

....

.... fKx 2,Ky 1... fKx 2,Ny 1 fKx 2,0... fKx 2,Ky Тогда для матрицы def G(Nx,Ny ) = (F(Nx,Ny ) )1 diag(F(Nx,Ny ) vec G)F(Nx,Ny ) выполняются равенства W(Lx,Ly ) (F) 0Kx Ly (Nx,Ny ), (5.3.1) vec = vec G 0Lx 1Ky 0Lx 1Ly W(Kx,Ky ) (F) 0Kx 1Ky 1 0Lx Ky (Nx,Ny ). (5.3.2) vec = vec G 0Kx 1Ly Доказательство.

В силу линейности равенств, достаточно доказать равенства (5.3.1) и (5.3.2) для F = F (1) (F (2) )T и = U (1) (V (2) )T. Действительно, по предложениям 1.1.1, 1.1.2 и 1.1.3 выполняется 0Kx,Ly G(Nx,Ny ) vec = 0Lx 1,Ky 0Lx 1,Ly = (F1 F1 )diag((FNy FNx )(F (2) F (1) )) Ny Nx V (2) V (1) (FNy FNx ) = 0Ly 1 0Lx V (2) V (1) F1 diag(FNy F (2) )FNy F1 diag(FNx F (1) )FNx = = Ny Nx 0Ly 1 0Lx X(Ly ) (FNy )V (2) X(Lx ) (FNx )V (1) = = W(Lx,Ly ) (F)V (1) (V (2) )T = vec, что и требовалось доказать.

5.3.3. Структура и краткое описание комплекса программ Комплекс программ реализован, в основном на С++. Объем комплекса программ составляет более 90 классов, 28000 строк, более 100 файлов.

Комплекс программ состоит из двух независимых частей.

• Набор программ, основанных на ядре 2DSSA, реализующего алгоритмы из раздела 5.3.1.

• Модуль пакета RSSA для двумерного АСС, реализующий методы из раз дела 5.3.2 и набор программ на языке R.

Пакет RSSA, основным разработчиком которого является А.И. Коробей ников [81], реализует набор функций для АСС-разложения рядов и массивов в среде для математических вычислений R [53]. Модуль для двумерного АСС, разработанный соискателем, позволяет осуществлять быстрое разложение мас сивов с помощью БПФ. Многие результаты раздела 5.1 получены с помощью программ на языке R, использующих данный модуль.

Ядро 2DSSA представляет собой кроссплатформенную библиотеку на язы ке С++, предоставляющую программный объектно-ориентированный интер фейс для АСС-разложения массивов с помощью вычисления ковариационной матрицы. На ядре основаны программы:

интерактивная программа для ОС Windows, предоставляющая с • 2dssa помощью оконного интерфейса наглядное представление исходного мас сива, собственных и факторных массивов, и т.д. Программа позволяет полностью контролировать процесс вычислений. Программа может за гружать и сохранять из файлов текстовых форматов и наиболее распро страненных форматов цифровых изображений.

консольная программа для тестирования очистки сигна • modelling ла от шума методом АСС. Программа позволяет проводить численные эксперименты по нахождению оптимального количества компонент, оп тимальных размеров окна, сравнивать качество восстановления для раз личных областей, моделировать шумы разного типа.

консольная программа для задач обработки текстурных изоб • texture ражений. Программа состоит из нескольких подпрограмм для задач клас сификации и сегментации с помощью методов, описанных в разделе 5.2.2.

Результаты раздела 5.2.2 получены с помощью данной программы.

Заключение В диссертационной работе исследовались проблемы разделимости времен ных и пространственных данных в методе АСС и описания структуры моделей рядов и массивов конечного ранга, порождаемых разложением в методе АСС.

Были предложены новые математические методы решения данных проблем.

Были построены алгоритмы для решения задач обработки двумерных данных, свойства моделей и алгоритмов были исследованы методами статистического моделирования. Алгоритмы были эффективно реализованы в виде комплекса программ. Рассмотрим подробнее основные выводы работы.

• Было установлено, что в рамках алгебраической теории ганкелевых мат риц становится возможным систематизировать и установить точные свя зи между типами рядов: конечного ранга, реверсивными, продолжимы ми, и т.д. (предложения 3.1.3, 3.1.4, 3.1.6 и 3.1.7;

теорема 3.1.1), а также систематизировать важные для практики свойства, такие как свойства побочных корней ЛРФ прогноза (см. теорему 3.3.1, предложение 3.3.3 и раздел 3.3.3).

• Был предложен новый критерий полуразделимости (теорема 3.2.1), кото рый позволил найти все случаи полуразделимости и разделимости про извольных рядов (теоремы 3.2.3 и 3.2.4) и разделимости массивов конеч ного ранга (теорема 4.3.3). Данный критерий может быть полезен для дальнейшего исследования разделимости в АСС, например, более слож ного, практически и теоретически значимого случая асимптотической разделимости.

• Был предложен способ нахождения количества ганкелевых матриц дан ного ранга и определенной структуры (следствие 3.4.3, предложения 3.4. и 3.4.5), которое необходимо для вычисления вероятностей линейной классификации с помощью ганкелевых матриц в варианте метода АСС над конечным полем. Полученные результаты в дальнейшем могут быть использованы для нахождения границ доверительных областей в соот ветствующих статистических тестах.

• Было показано, что бесконечные массивы конечного ранга соответству ют массивам с конечной линейной сложностью (линейным рекуррент ным массивам, см. предложение 4.1.1);

был найден общий вид непре рывных функций, имеющих конечное АСС разложение в непрерывном варианте метода АСС (теорема 4.4.2).

• В рамках исследования была получена новая оценка линейной слож ности (ранга) полиномиального массива по его коэффициентам (теоре ма 4.1.2), которая в некоторых случаях лучше известной оценки (предло жение 2.2.6). Однако и полученную оценку, возможно, можно улучшить.

Оценка ранга массива важна для предварительного этапа математиче ского моделирования в рамках метода АСС.

• Был исследован вопрос о нахождении области допустимых размеров ок на для массивов конечного ранга. На общий случай массивов конечного ранга расширены оценки множества допустимых размеров окна с помо щью методов базисов Грёбнера (теорема 4.2.1). Знание минимально до пустимых размеров окна необходимо для выбора достаточных размеров окна в практических задачах.

• Была исследована зависимость разделимости сигнала и шума от разме ров окна для различных массивов конечного ранга. Было продемонстри ровано, что многие эффекты, известные и доказанные для одномерного случая, имеют место и для двумерного АСС, и в некоторых случаях усиливаются. Было показано, что двумерный АСС в задачах восстанов ления зашумленного сигнала во многих случаях предпочтительнее обыч ного сингулярного разложения, которое используется во многих методах обработки данных (см. раздел 5.1).

• Метод АСС был применен к задаче фильтрации цифровых моделей ре льефа. Были разработаны методы обработки текстурных изображений на основе введения расстояния, что упрощает процедуру классификации по сравнению с известными методами, основанными на АСС (см. раз дел 5.2.2). На тестовых изображениях и базах данных изображений были продемонстрированы эффективность разработанных методов и преиму щества использования больших размеров окна.

• Были разработаны и реализованы в виде комплекса программ алгорит мы для двумерного разложения в методе АСС и для решения различных задач обработки данных (см. раздел 5.3).

Таким образом, в диссертации были отражены теоретические и методологиче ские вопросы математического моделирования и разработки устойчивых алго ритмов, реализованных в виде комплекса программ, пригодных для решения научных и практических задач.

Литература 1. Александров Ф. И. Выделение аддитивных компонент временного ряда при пакетной обработке методом “Гусеница”-SSA // Вестник СПбГУ, Cер.

1. Математика. Механика. Астрономия. 2006. № 2. С. 71–74.

2. Алексеева Н. Комбинаторный анализ двух форм скрытой периодичности категориальных последовательностей // Вестник Санкт-Петербургского университета. Серия 1: Математика, механика, астрономия. 2007. № 3.

С. 55–64.

3. Андерсон Т. Статистический анализ временных рядов. М.: Мир, 1976.

С. 757.

4. Антоновский М., Бухштабер В., Векслер Л. Применение многомерно го статистического анализа для обнаружения структурных изменений во временных рядах данных экологических наблюдений // Проблемы экологического мониторинга и моделирования экосистем. 1993. Т. 15.

С. 193–212.

5. Бардин М. Изменчивость температуры воздуха над западными террито риями России и сопредельными странами в XX веке // Метеорология и гидрология. 2002. № 8. С. 5–23.

6. Белонин М. Д., Татаринов И. В., Калинин О. М. и др. Факторный анализ в нефтяной геологии. М.: ВИЭМС, 1971. С. 56.

7. Бриллинджер Д. Временные ряды. Обработка данных и теория. М.: Мир, 1980. С. 536.

8. Бухштабер В. М. Многомерные развертки временных рядов. Теоретиче ские основы и алгоритмы // Обозрение прикл. промышл. матем., сер.

Вероятн. и статист. 1997. Т. 4, № 4. С. 629–645.

9. Ван дер Варден Б. Л. Алгебра. М.: Наука, 1976. С. 623.

10. Владимиров В. С. Уравнения математической физики. М.: Наука, 1981.

С. 512.

11. Воеводин В., Кузнецов Ю. Матрицы и вычисления. М.:Наука, 1984.

С. 320.

12. Гантмахер Ф. Р. Теория матриц, 5-е изд. ФИЗМАТЛИТ, 2004. С. 576.

13. Главные компоненты временных рядов: метод Гусеница, Под ред.

Д. Л. Данилов, А. А. Жиглявский. СПб.: Пресском, 1997. С. 308.

14. Глухов М. М., Елизаров В. П., Нечаев А. А. Алгебра, Учебник. В 2-х т.

Т. II. М.: Гелиос АРВ, 2003.

15. Голуб Д., Лоан Ч. В. Матричные вычисления. М.: Мир, 1999. С. 548.

16. Голяндина Н. Э. Метод Гусеница -SSA: анализ временных рядов: Учеб.

пособие. СПб.: ВВМ, 2003. С. 85.

17. Голяндина Н. Э. Метод Гусеница -SSA: прогноз временных рядов: Учеб.

пособие. СПб.: ВВМ, 2004. С. 53.

18. Голяндина Н. Э., Осипов Е. В. Метод “Гусеница”-SSA для анализа вре менных рядов с пропусками // Мат. модели. Теория и приложения. 2005.

6. С. 50–61.

19. Голяндина Н. Э., Усевич К. Д. Метод 2D-SSA для анализа двумерных по лей // Труды VII Межд. конф Идентификация систем и задачи управ ления SICPRO’08. Москва, 28–31 января 2008. 2008. С. 1657–1727.

20. Голяндина Н. Э., Флоринский И. В., Усевич К. Д. Анализ сингулярно го спектра для фильтрации цифровых моделей рельефа // Геодезия и картография. 2008. № 5. С. 21–28.

21. Дюк В. Поиск сложных непериодических шаблонов в последовательно стях числе и символов методами локальной геометрии // Тр. СПИИ РАН.

2002. Т. 2, № 1. С. 263–268.

22. Ефимов В. М., Галактионов Ю. К. О возможности прогнозирования цик лических изменений численности млекопитающих // Журн. общ. биоло гии. 1983. Т. 44, № 3. С. 343–352.

23. Иохвидов И. С. Ганкелевы и тeплицевы матрицы и формы. Наука, М., 1974. С. 263.

24. Истомин И. А., Котляров О. Л., Лоскутов А. Ю. К проблеме обработки временных рядов: расширение возможностей метода локальной аппрок симации посредством сингулярного спектрального анализа // Теоретиче ская и математическая физика. 2005. Т. 142, № 1. С. 148–159.

25. Кислицин М. М. Многомерная статистика временных рядов наблюдений в авиационной эргономике // Вопросы кибернетики.Биотехнические си стемы в авиационной эргономике. 1978. № 51. С. 117–126.

26. Кокс Д., Литтл Д., О’Ши Д. Идеалы, многообразия и алгоритмы. М.:

Мир, 2000. С. 687.

27. Кузьмин А. С., Куракин В. Л., Нечаев А. А. Псевдослучайные и полили нейные последовательности // Труды по дискретной математике / Под ред. П. ред. В.Н. Сачкова и др. М.: ТВП, 1997. Т. 1. С. 139–202.

28. Куракин В. Л. Линейная сложность полилинейных последовательно стей // Дискретная математика. 2001. Т. 13. С. 4–55.

29. Куракин В. Л. Биномиальная сложность полилинейных последовательно стей // Труды по дискр. матем. М.: ФИЗМАТЛИТ, 2002. Т. 6. С. 82–138.

30. Лидл Р., Нидеррайтер Г. Конечные поля: в 2-х т. Государственное издание физико-математической литературы, 1988. С. 820.

31. Линник Ю. Метод наименьших квадратов и основы теории обработки наблюдений. М.: Физматлит, 1958. С. 336.

32. Марпл-мл. С. Л. Цифровой спектральный анализ и его приложения. М.:

МИР, 1990. С. 584.

33. Муттер В. М. Основы помехоустойчивой телепередачи информации. Л.:

Энергоатомиздат, 1990. С. 288.

34. Некруткин В. Аппроксимирующие пространства и продолжения времен ных рядов // Статистические модели с приложениями в эконометрике и смежных областях / Под ред. П. ред. С.М.Ермакова и Ю.Н.Каштанова.

СПб.: изд-во НИИХ СПбГУ, 1999. С. 2–32.

35. Некруткин В. В. Разложения временных рядов // Главные компоненты временных рядов: метод Гусеница / Под ред. А. А. Ж. Под ред. Д.

Л. Данилова. СПб: Пресском, 1997. С. 194–227.

36. Прасолов В. В. Многочлены. М.: МЦНМО, 2001. С. 336.

37. Прэтт У. Цифровая обработка изображений, Под ред. пер. с англ.

под ред. Д. С. Лебедева. М. : Мир, 1982. С. 310.

38. Усевич К. Д. Разложение функций в двумерном варианте метода Гусени ца -SSA и связанные с ним системы уравнений в частных производных // Вестник СПбГУ, Серия 10: Прикладная математика, информатика, про цессы управления. 2009. № 3. С. 152–161.

39. Ade F. Characterization of Textures by ‘Eigenlters’ // Signal Processing.

1983. Vol. 5. Pp. 451–457.

40. Alexeyeva N., Alexeyev A., Gracheva P. et al. Symptom and syndrome analysis of categorial series, logical principles and forms of logic // Proc. of the 2010 3rd Int. Сonf. on BioMed. Engineering and Informatics.

Vol. 6. 2010. Pp. 2603–2606.

41. Althaler J., Dur A. Finite linear recurring sequences and homogeneous ide als // Applicable Algebra in Engineering, Communication and Computing.

1996. Vol. 7. Pp. 377–390.

42. Axmon J., Hansson M., Sornmo L. Partial forward-backward averaging for enhanced frequency estimation of real X-texture modes // IEEE Transactions on Signal Processing. 2005. Vol. 53, no. 7. Pp. 2550–2562.

43. Badeau R., David B., Richard G. High-Resolution Spectral Analysis of Mix tures of Complex Exponentials Modulated by Polynomials // IEEE Transac tions on Signal Processing. 2006. Vol. 54, no. 4. Pp. 1341–1350.

44. Bezerra L. H., Bazan F. S. V. Eigenvalue Locations of Generalized Companion Predictor Matrices // SIAM Journal on Matrix Analysis and Applications.

1998. Vol. 19, no. 4. Pp. 886–897.

45. Biondi F., Isaacs C., Hughes M. et al. The Near-1600 Dry/wet Knock-out:

Linking Terrestrial and Near-shore Ecosystems // Proc. 24th Annual Cli mate Diagnostics and Prediction Workshop, / US department of Commerce, Springeld, VA. 76–79.

46. Broomhead D. S., King G. P. Extracting qualitative dynamics from experi mental data // Physica. 1986. Vol. 20. Pp. 217–236.

47. Buchstaber V. M. Time Series Analysis and Grassmanians // Applied prob lems of Radon transform / Ed. by S. Gindikin. 1994. Vol. 162 of American Mathematical Society Translations. Pp. 1–18.

48. Chateld C. The Analysis of Time Series: An Introduction. 2nd edition.

Chapman&Hall, London, 1980. P. 268.

49. Cheveigne A. D., Simon J. Z. Denoising based on time-shift PCA // Journal of Neuroscience Methods. 2007. Vol. 165. Pp. 297–305.

50. Chi-Jie L., Du-Ming T. Automatic defect inspection for LCDs using singular value decomposition // The International Journal of Advanced Manufactur ing Technology. 2005. Vol. 25. Pp. 53–61.

51. Colebrook J. M. Continuous plankton records zooplankton and environ ment, northeast Atlantic and North Sea // Oceanol. Acta, 1948-1975. 1978.

Vol. 1. Pp. 9–23.

52. Columbia-Utrecht Reectance and Texture Database (CUReT). 1999.

URL: http://www1.cs.columbia.edu/CAVE/software/curet/ (дата обра щения: 08.09.2010).

53. Dalgaard P. Introductory Statistics with R. 2nd edition. Springer, 2008.

P. 380.

54. Dana K. J., van Ginneken B., Nayar S. K., Koenderink J. J. Reectance and Texture of Real World Surfaces // ACM Transactions on Graphics. 1999.

Vol. 18, no. 1. Pp. 1–34.

55. Daykin D. E. Distribution of bordered persymmetric matrices in a nite eld // J. Reine und Angew. Math. (Crelles Journal). 1960. Vol. 203.

Pp. 47–45.

56. Elkies N. D. On nite sequences satisfying linear recursions // New York Journal of Mathematics. 2002. Vol. 8. P. 85–97.

57. Feldmann S., Heinig G. Parametrization of minimal rank block Hankel ma trix extensions and minimal partial realizations // Integral Equations and Operator Theory. 1999. Vol. 33. Pp. 151–171.

58. Findley D. F. On the Early History of the Singular Value Decomposition // SIAM Review. 1993. Vol. 35, no. 4. Pp. 551–566.

59. Florinsky I. Derivation of topographic variables from a digital elevation model given by a spheroidal trapezoidal grid // International Journal of Geograph ical Information Science. 1998. Vol. 12. Pp. 829–852.

60. Florinsky I. Errors of signal processing in digital terrain modelling // Inter national Journal of Geographical Information Science. 2002. Vol. 16, no. 5.

Pp. 475–501.

61. Ghil M., Allen R., Dettinger M. et al. Advanced spectral methods for climatic time series // Rev. Geophys. 2002. Vol. 40, no. 1. Pp. 1–41.

62. Ghil M., Vautard R. Interdecadal oscillations and the warming trend in global temperature time series // Nature. 1991. Vol. 350. Pp. 324–327.

63. Golub G., Reinsch C. Singular value decomposition and least squares solu tions // Numerical Mathematics. 1970. Vol. 14. Pp. 403–420.

64. Golyandina N. On the choice of parameters in Singular Spectrum Analysis and related subspace-based methods // Statistics and Its Interface. 2010.

Vol. 3, no. 3. Pp. 259–279.

65. Golyandina N., Florinsky I., Usevich K. Filtering of Digital Terrain Models by Two Dimensional Singular Spectrum Analysis // International Journal of Ecology & Development. 2007. Vol. 8, no. F07. Pp. 81–94.

66. Golyandina N., Nekrutkin V., Zhigljavsky A. A. Analysis of Time Series Struc ture: SSA and Related Techniques. Chapman and Hall/CRC, 2001. P. 320.

67. Golyandina N., Osipov E. The “Caterpillar”-SSA method for analysis of time series with missing values // J. Stat. Plan. Infer. 2007. Vol. 137, no. 8.

Pp. 2642–2653.

68. Golyandina N., Usevich K. An Algebraic View on Finite Rank in 2D-SSA // Proceedings of the 6th St.Petersburg Workshop on Simulation, June-July 2009 / Ed. by S. Ermakov, V. Melas, A. Pepelyshev. 2009. Pp. 308–313.

69. Golyandina N., Vlassieva E. First-order SSA-errors for long time series: model examples of simple noisy signals // Proceedings of the 6th St.Petersburg Workshop on Simulation Vol.1, June 28-July 4, 2009, St. Petersburg. 2009.

Pp. 314–319.

70. Golyandina N. E., Usevich K. D. 2D-extension of Singular Spectrum Analysis:

algorithm and elements of theory // Matrix methods: theory, algorithms and applications / Ed. by V. Olshevsky, E. Tyrtyshnikov. World Scientic, 2010.

Pp. 449–473.

71. Govil N. K., Rahman Q. I. On the Enestrm-Kakeya theorem. // Tohoku o Mathematical Journal. 1968. Vol. 20, no. 2. Pp. 126–136.

72. Hakopian H., Tonoyan M. Partial dierential analogs of ordinary dierential equations and systems // New York Journal of Mathematics. 2004. Vol. 10.

Pp. 89–116.

73. Handbook of texture analysis, Ed. by M. Mirmehdi, X. Xie, J. Suri. Imperial College Press, 2008. P. 423.

74. Heinig G., Rost K. Algebraic methods for Toeplitz-like matrices and opera tors. Akademie Verlag, Berlin, 1984. P. 212.

75. Holloway D. M., Harrison L. G., Kosman D. et al. Analysis of Pattern Pre cision Shows That Drosophila Segmentation Develops Substantial Indepen dence From Gradients of Maternal Gene Products // Development Dynamics.

2006. Vol. 235. Pp. 2949–2960.

76. Hsieh W. W., Hamilton K. Nonlinear singular spectrum analysis of the trop ical stratospheric wind // Quarterly Journal of the Royal Meteorological So ciety. 2003. Vol. 129. Pp. 2367–2382.

77. Jevrejeva S., Moore J. C. Singular Spectrum Analysis of Baltic Sea ice condi tions and large-scale atmospheric patterns since 1708 // Geophy. Res. Lett.

2001. Vol. 28, no. 23. Pp. 4503–4506.

78. Johnson S. G., Frigo M. A modied split-radix FFT with fewer arithmetic op erations // IEEE Trans. Signal Processing. 2007. Vol. 55, no. 1. Pp. 111–119.

79. Kehrein A., Kreuzer M. Computing Border Bases // Journal of Pure and Applied Algebra. 2006. May. Vol. 205, no. 2. Pp. 279–295.

80. Konstantinides K., Natarajan B., Yovanof G. S. Noise estimation and lter ing using block-based singular value decomposition // IEEE Transactions on Image Processing. 1997. Vol. 6. Pp. 479–483.

81. Korobeynikov A. Computation- and space-ecient implementation of SSA // Statistics and Its Interface. 2010. Vol. 3, no. 3. Pp. 357–368.

82. Krim H., Forster P., Proakis J. G. Operator approach to performance analysis of root-MUSIC and root-min-norm // IEEE Transactions on Signal Process ing. 1992. Vol. 40, no. 7. Pp. 1687–1696.

83. Kumaresan R. On the Zeros of the Linear Prediction-Error Filter for De terministic Signals // IEEE Transactions on Acoustics, Speech, and Signal Processing. 1983. Vol. 31, no. 1. Pp. 217–220.

84. Kumaresan R., Tufts D. W. Data-adaptive principal component signal pro cessing // Proc. of IEEE Conference On Decision and Control. Albuquerque, 1980. Pp. 949–954.

85. Kumaresan R., Tufts D. W. Estimating the Parameters of Exponentially Damped Sinusoids and Pole-Zero Modelling in Noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. 1982. Vol. 30, no. 6. Pp. 833–840.

86. Kumaresan R., Tufts D. W. Estimating the Angles of Arrival of Multiple Plane Waves // IEEE Transactions on Aerospace and Electronic Systems.

1983. Vol. AES-19, no. 1. Pp. 134–139.

87. Laub A. J. Matrix Analysis for Scientists and Engineers. SIAM, 2004.

88. Laws K. I. Textured image segmentation: Report 940 / Image Processing Institute, University of southern California. 1980. P. 190.

89. Li F., Liu H., Vaccaro R. J. Performance analysis for DOA estimation algo rithms: unication, simplication, and observations // IEEE Transactions on Aerospace and Electronic Systems. 1993. Vol. 29, no. 4. Pp. 1170–1184.

90. Loan C. F. V., Pitsianis N. P. Approximation with Kronecker products // Lin ear Algebra for Large Scale and Real Time Applications / Ed. by M. S. Moo nen, e. G. H. Golub. Kluwer Publications, 1993. Pp. 293–314.

91. Magnus J. R., Neudecker H. Matrix Dierential Calculus with Applications to Statistics and Econometrics. John Wiley & Sons, 2004. P. 450.

92. Mahecha M., Reichstein M., Lange H. et al. Characterizing ecosystem-atmo sphere interactions from short to interannual timescales // Biogeosciences Discuss. 2007. Vol. 4. Pp. 1405–1435.

93. Mamou J., Feleppa E. J. Singular spectrum analysis applied to ultrasonic detectionand imaging of brachytherapy seeds // J. Acoust. Soc. Am. 2007.

Vol. 121, no. 3. Pp. 1790–1801.

94. Markovic D., Koch M. Characteristic scales, temporal variability modes and simulation of monthly Elbe River ow time series at ungauged locations // Physics and Chemistry of the Earth. 2006. Vol. 31. Pp. 1262–1273.

95. Mart inez-Finkelshtein A., McLaughlin K. T.-R., Sa E. B. Asymptotics of orthogonal polynomials with respect to an analytic weight with algebraic singularities on the circle // International Mathematics Research Notices.

2006. Vol. 2006, no. 91426. Pp. 1–43.

96. Mart inez-Finkelshtein A., McLaughlin K. T.-R., Sa E. B. Szeg Orthogonal o Polynomials with Respect to an Analytic Weight: Canonical Representation and Strong Asymptotics // Constr. Approx. 2006. Vol. 24, no. 3. Pp. 319–363.

97. McElroy T., Sutclie A. An iterated parametric approach to nonstationary signal extraction // Computational Statistics & Data Analysis. 2006. Vol. 50.

Pp. 2206–2231.

98. Monadjemi A. Towards Ecient Texture Classication and Abnormality De tection: Ph. D. thesis / Department of Computer Science, University of Bris tol. 2004. P. 177.

99. Monadjemi A., Mirmehdi M., Thomas B. Restructured eigenlter matching for novelty detection in random texture // In Proceedings of the 15th British Machine Vision Conference. 2004. Pp. 637–646.

100. Mourrain B. A New Criterion for Normal Form Algorithms // Proceedings of the 13th International Symposium on Applied Algebra, Algebraic Algorithms and Error-Correcting Codes. AAECC-13. 1999. Pp. 430–443.

101. Nekrutkin V. Perturbation expansions of signal subspaces for long signals // Statistics and Its Interface. 2010. Vol. 3, no. 3. Pp. 297–319.

102. Paegle J. N., Byerle L., Mo K. Intraseasonal Modulation of South American Summer // Monthly Weather Review. 2000. March. Vol. 128. Pp. 837–850.

103. Pakula L. Asymptotic Zero Distribution of Orthogonal Polynomials in Sinu soidal Frequency Estimation // IEEE Transactions on Information Theory.

1987. Vol. 33, no. 4. Pp. 569–576.

104. Pan V. Structured matrices and polynomials: unied superfast algorithms.

Birkhuser Boston, 2001. P. 278.

a 105. Patel D., Davies E., Hannah I. The use of convolution-operators for detecting contaminants in food images // Pattern Recognition. 1996. Vol. 29, no. 6.

Pp. 1019–1029.

106. Planas C. The Analysis of Seasonality in Economic Statistics: A Survey of Re cent Developments: Tech. rep.: EUROSTAT working group document, 1997.

107. Randen T., Husoy J. H. Filtering for Texture Classication: A Comparative Study // IEEE Transactions on Pattern Analysis and Machine Intelligence.

1999. Vol. 21, no. 4. Pp. 291–310.

108. Rezek I., Roberts S. Stochastic Complexity Measures for Physiological Sig nal Analysis // IEEE Transactions on BME. 1998. Sept. Vol. 25, no. 9.

Pp. 1186–1191.

109. Rodriguez-Aragon L. J., Zhigljavsky A. Singular spectrum analysis for image processing // Statistics and Its Interface. 2010. Vol. 3, no. 3. Pp. 419–426.

110. Rouquette S., Najim M. Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods // IEEE Transactions on Signal Pro cessing. 2001. Vol. 49, no. 1. Pp. 237–245.

111. Roy R., Kailath T. ESPRIT: estimation of signal parameters via rotational invariance techniques // IEEE Trans. Acoust. 1989. Vol. 37. Pp. 984–995.

112. Serita A., Hattori K., Yoshino C., M. Hayakawa N. I. Principal component analysis and singular spectrum analysis of ULF geomagnetic data associat ed with earthquakes // Natural Hazards and Earth System Sciences. 2005.

Vol. 5. Pp. 685–689.

113. Simon B. Orthogonal Polynomials on the Unit Circle, Part 1: Classical The ory. AMS Colloquium Series, AMS, Providence, RI, 2005. P. 466.

114. SRTM3 (Shuttle Radar Topographic Mission) Version 2. 2003. URL:

http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/ (дата обращения:

08.10.2011).

115. Stepanov D., Golyandina N. SSA-based approaches to analysis and forecast of multidimensional time series // Proceedings of the 5th St.Petersburg Work shop on Simulation, St. Petersburg State University, St. Petersburg. 2005.

Pp. 293–298.

116. Stetter H. J. Numerical polynomial algebra. SIAM, Philadelphia, 2004. P. 472.

117. Stoica P., Moses R. L. Introduction to spectral analysis. Prentice Hall, 1997.

P. 319.

118. Strang G. Introduction to Linear Algebra. 2003. P. 568.

119. Szabados J. On some problems connected with orthogonal polynomials // Acta Mathematica Academiae Scientarium Hungaricae. 1979. Vol. 33, no.

1–2. Pp. 197–210.

120. Thomakos D., Wanga T., Wille L. Modeling daily realized futures volatility with singular spectrum analysis // Physica A: Statistical Mechanics and its Applications. 2002. Vol. 312, no. 3–4. Pp. 505–519.

121. Unser M., Ade F. Feature extraction and decision procedure for automated inspection of textured materials // Pattern Recognition Letters. 1984.

March. Vol. 2. Pp. 181–195.

122. Usevich K. On signal and extraneous roots in Singular Spectrum Analysis // Statistics and Its Interface. 2010. Vol. 3, no. 3. Pp. 281–295.

123. Varma M., Zisserman A. Texture Classication: Are Filter Banks Neces sary? // IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 18–20 June 2003. Proceedings. Vol. 2. 2003. Pp. 691–698.


124. Varma M., Zisserman A. A Statistical Approach to Texture Classication from Single Images // Int. J. Comput. Vision. 2005. Vol. 62, no. 1–2.

Pp. 61–81.

125. Vassiliev N., Pavlov D. Strong non-Noetherity of polynomial reduction // Записки научных семинаров ПОМИ. 2009. Т. 373. С. 73–76.

126. Wang F., Wang S., Dou H. Estimating two-dimensional frequency pairs by extended Esprit-type method // Proceedings, ICCT 2003. International Con ference on Communication Technology. Vol. 2. 2003. Pp. 1464–1467.

127. Yang H. H., Hua Y. On Rank of Block Hankel Matrix for 2-D Frequency Detection and Estimation // IEEE Transactions on Signal Processing. 1996.

Vol. 44, no. 4. Pp. 1046–1048.

128. Yu X.-L., Buckley K. M. Bias and variance of direction-of-arrival estimates from MUSIC, MIN-NORM, and FINE // IEEE Transactions on Signal Pro cessing. 1994. Vol. 42, no. 7. Pp. 1812–1816.

Приложение А Исходные коды комплекса программ А.1. Ядро комплекса программ 2D-SSA Файл “eldssabasiccalc.cpp”:

#include "comdefs.h" #include sstream #include memory.h #include string using namespace std;

#include math.h #include "fieldssabasic.h" #include "libs\clapack\clapack.h" #define CovTr(u1, v1, u2, v2) \ (XXt[(((v1) * Lx + (u1)) * ((v1) * Lx + (u1) + 1)) / 2 + \ (v2) * Lx + (u2)]) #define CovSq(u1, v1, u2, v2) \ (XXt[((v1) * Lx + (u1)) * L + (v2) * Lx + (u2)]) inline void SetCovValue ( real *XXt, int Lx, int L, bool triangular, int u1, int v1, int u2, int v2, real val ) { if (triangular) { CovTr(u1,v1,u2,v2) = val;

} else { CovSq(u1,v1,u2,v2) = CovSq(u2,v2,u1,v1) = val;

} } inline void AddCovValue ( real *XXt, int Lx, int L, bool triangular, int u1, int v1, int u2, int v2, int u1src, int v1src, int u2src, int v2src, real val ) { if (triangular) { CovTr(u1,v1,u2,v2) = val + CovTr(u1src,v1src,u2src,v2src);

} else { CovSq(u1,v1,u2,v2) = CovSq(u2,v2,u1,v1) = val + CovSq(u1src,v1src,u2src,v2src);

} } void FieldSSABasic::CalculateCovarMatrix( real *XXt, int Lx, int Ly, FieldNoAlloc &Fld, FieldNoAlloc *CentrComp, TaskInfoHandler *hand, bool triangular ) { /*** 2: Make covariations matrix ***/ int L = Lx * Ly, Nx = Fld.getNx(), Ny = Fld.getNy();

int Kx = Nx - Lx + 1, Ky = Ny - Ly + 1;

int K = Kx * Ky;

if (hand) hand-startNextTask(string("Make covariations matrix"), (L * (L + 1)) / 2);

{ int dy, y1, y2;

int dx, x1, x2;

int cnt = 0;

for (dy = 0;

dy Ly;

dy++) { /* Fill pivot (fully calculated) elements */ for (dx = 0;

dx Lx;

dx++) { SetCovValue(XXt, Lx, L, triangular, dx, dy, 0, 0, Fld.subFieldDotSubField(dx, dy, Fld, 0, 0, Kx, Ky) / K);

if (hand) hand-completedIterations(++cnt);

if (dy != 0 && dx != 0) { SetCovValue(XXt, Lx, L, triangular, 0, dy, dx, 0, Fld.subFieldDotSubField(0, dy, Fld, dx, 0, Kx, Ky) / K);

if (hand) hand-completedIterations(++cnt);

} } /* Fill first strip (calculated by shift by x) */ for (x1 = 1;

x1 Lx;

x1++) { int x_up_bound = (dy == 0 ? x1 : Lx - 1);

for (x2 = 1;

x2 = x_up_bound;

x2++) { AddCovValue(XXt, Lx, L, triangular, x1, dy, x2, 0, x1 - 1, dy, x2 - 1, 0, -(Fld.subFieldDotSubField(x1-1, dy, Fld, x2-1, 0, 1, Ky) Fld.subFieldDotSubField(x1-1+Kx, dy, Fld, x2-1+Kx, 0, 1, Ky)) / K);

if (hand) hand-completedIterations(++cnt);

} } /* Fill other strips based on first */ for (y2 = 1;

y2 Ly - dy;

y2++) { y1 = y2 + dy;

for (x1 = 0;

x1 Lx;

x1++) { int x_up_bound = (dy == 0 ? x1 : Lx - 1);

for (x2 = 0;

x2 = x_up_bound;

x2++) { AddCovValue(XXt, Lx, L, triangular, x1, y1, x2, y2, x1, y1 - 1, x2, y2 - 1, -(Fld.subFieldDotSubField(x1, y1 - 1, Fld, x2, y2 - 1, Kx, 1) Fld.subFieldDotSubField(x1, y1 - 1 + Ky, Fld, x2, y2 - 1 + Ky, Kx, 1)) / K);

if (hand) hand-completedIterations(++cnt);

} } } } if (CentrComp != NULL) { int i,j;

for (i = 0;

i L;

i++) { double mean_i = CentrComp-getValue(i%Lx, i/Lx), mean_j;

for (j = 0;

j = i;

j++) { mean_j = CentrComp-getValue(j%Lx, j/Lx);

AddCovValue(XXt, Lx, L, triangular, i%Lx, i/Lx, j%Lx, j/Lx, i%Lx, i/Lx, j%Lx, j/Lx, - mean_i * mean_j);

} } } } } void FieldSSABasic::SVDTransform( int lx, int ly, bool centring, SVD_METHOD method, int calc_eigen, TaskInfoHandler *hand ) { if (method != TRLAN && calc_eigen != 0) { calc_eigen = lx * ly;

} else { if (calc_eigen 0 || calc_eigen lx * ly) throw FieldSSAIllegalIndex(calc_eigen, lx * ly);

} if ((SvdMethod != method) || (Lx != lx) || (Ly != ly) || (Centring != centring)) { SetParametersAndCreateTriplesArrays(lx, ly, method);

} Centring = centring;

if (calc_eigen CalculatedEigens) return;

// We are done EigenFields-resizeArray(calc_eigen + 1);

real *XXt = NULL;

double *Ev = NULL;

try { if (hand) hand-setTotalTasks(3);

/* Calculate mean vector if is needed */ if (Centring) { CalculateCentrComp(&EigenFields-getItem(0), Lx, Ly, Fld);

} if (method == CLAPACK || method == JACOBI || method == TRLAN) { if (method == TRLAN) { XXt = new real[L * L];

CalculateCovarMatrix(XXt, Lx, Ly, Fld, (Centring ? &UncheckedGetEigenField(1) :

(FieldNoAlloc *)NULL), hand, false);

} else { XXt = new real[(L * (L + 1)) / 2];

CalculateCovarMatrix(XXt, Lx, Ly, Fld, (Centring ? &UncheckedGetEigenField(1) :

(FieldNoAlloc *)NULL), hand, true);

} Ev = new double[L+1];

Ev[0] = 0;

if (method != CLAPACK) { JacobiSVDCalculator calc;

calc.calculate(XXt, EigenFields-getDataPtr(Centring ? 1 : 0), L, L, hand, calc_eigen 0);

} else { ClapackSVDCalculator calc;

calc.calculate(XXt, EigenFields-getDataPtr(Centring ? 1 : 0), &Ev[(GetCentring() ? 1 : 0)], L, L, hand, calc_eigen 0);

} /*** 3.1: Convert eigen vectors to field ***/ if (hand) hand-startNextTask(string("Convert eigen" "vectors to field"), L);

real *PLambda = Ev;

LambdaSum = 0;

for (int triple = 1;

triple = L;

triple++, PLambda ++) { /* Store eigen value */ Pairs[triple - 1]-Lambda = *PLambda;

Pairs[triple - 1]-LambdaSqrt = sqrt(*PLambda);

LambdaSum += *PLambda;

} if (calc_eigen 0 ) { for (int triple = 1;

triple = L;

triple++, PLambda ++) { FieldNoAlloc &CurEigen = UncheckedGetEigenField(triple);

CurEigen.recalculateMaxMinValues();

if (hand) hand-completedIterations(triple);

} } } else { hand-completedIterations(L);

} delete[] XXt;


delete Ev;

Ev = NULL;

XXt = NULL;

} catch (FieldSSAInterruptedException e) { &e;

if (XXt != NULL) { delete XXt;

XXt = NULL;

} if (Ev != NULL) { delete [] Ev;

Ev = NULL;

} throw;

} /*** 5: Calculate percentage ***/ try { real sum = 0.0;

for (int triple = 1;

triple = L;

triple++) { Pairs[triple - 1]-Percent = (Pairs[triple - 1]-Lambda / LambdaSum) * 100.0;

sum += Pairs[triple - 1]-Percent;

Pairs[triple - 1]-CumPercent = sum;

} } catch (FieldSSAInterruptedException e) { &e;

throw;

} CalculatedEigens = calc_eigen;

if (State STATE_TRANSFORMED) { State = STATE_TRANSFORMED;

} } А.2. Модуль пакета RSSA для двумерного АСС Файл “hbhankel.cpp”:

/* * R package for Singular Spectrum Analysis * Copyright (c) 2009-10 Anton Korobeynikov * asl@math.spbu.ru * Copyright (c) 2009-10 Konstantin Usevich * usevich.k.d@gmail.com * * This program is free software;

you can redistribute it * and/or modify it under the terms of the GNU General Public * License as published by the Free Software Foundation;

* either version 2 of the License, or (at your option) * any later version.

* * This program is distributed in the hope that it will be * useful, but WITHOUT ANY WARRANTY;

without even the implied * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR * PURPOSE. See the GNU General Public License for more details.

* * You should have received a copy of the GNU General Public * License along with this program;

if not, write to the * Free Software Foundation, Inc., 675 Mass Ave, Cambridge, * MA 02139, USA.

*/ #include R.h #include Rinternals.h #include complex.h #include fftw3.h #include "extmat.h" typedef struct { fftw_complex * circ_freq;

fftw_plan r2c_plan;

fftw_plan c2r_plan;

struct {R_len_t x;

R_len_t y;

} window;

struct {R_len_t x;

R_len_t y;

} length;

} hbhankel_matrix;

static unsigned hbhankel_nrow(const void *matrix) { const hbhankel_matrix *h = matrix;

return h-window.x * h-window.y;

} static unsigned hbhankel_ncol(const void *matrix) { const hbhankel_matrix *h = matrix;

return (h-length.x - h-window.x + 1)* (h-length.y - h-window.y + 1);

} static void free_circulant(hbhankel_matrix *h) { fftw_free(h-circ_freq);

fftw_destroy_plan(h-r2c_plan);

fftw_destroy_plan(h-c2r_plan);

} static void initialize_circulant(hbhankel_matrix *h, const double *F, R_len_t Nx, R_len_t Ny, R_len_t Lx, R_len_t Ly) { R_len_t Kx = Nx - Lx + 1, Ky = Ny - Ly + 1, i, j;

fftw_complex *ocirc;

fftw_plan p1, p2;

double *circ;

/* Allocate needed memory */ circ = (double*) fftw_malloc(Nx * Ny * sizeof(double));

ocirc = (fftw_complex*) fftw_malloc(Ny * (Nx/2 + 1) * sizeof(fftw_complex));

/* Estimate the best plans for given input length, note, that input data is stored in column-major mode, that’s why we’re passing dimensions in *reverse* order */ p1 = fftw_plan_dft_r2c_2d(Ny, Nx, circ, ocirc, FFTW_ESTIMATE);

p2 = fftw_plan_dft_c2r_2d(Ny, Nx, ocirc, circ, FFTW_ESTIMATE);

/* Fill input buffer */ /* TF - cbind(F[,Ky:Ny],F[,1:(Ky-1)]);

TF - rbind(TF[Kx:Nx,],TF[1:(Kx-1),]);

*/ for (j = 0;

j Ny;

++j) for (i = 0;

i Nx;

++i) /* This is pretty ad-hoc solution and needs to be fixed in the future */ circ[i + Nx*j] = F[(i+Kx-1) % Nx + Nx*((j+Ky-1) % Ny)];

/* Run the plan on input data */ fftw_execute(p1);

/* Cleanup and return */ fftw_free(circ);

h-circ_freq = ocirc;

h-r2c_plan = p1;

h-c2r_plan = p2;

h-window.x = Lx;

h-window.y = Ly;

h-length.x = Nx;

h-length.y = Ny;

} static void hbhankel_matmul(double* out, const double* v, const void* matrix) { const hbhankel_matrix *h = matrix;

R_len_t Nx = h-length.x, Ny = h-length.y;

R_len_t Lx = h-window.x, Ly = h-window.y;

R_len_t Kx = Nx - Lx + 1, Ky = Ny - Ly + 1, i, j;

double *circ;

fftw_complex *ocirc;

/* Allocate needed memory */ circ = (double*) fftw_malloc(Nx * Ny * sizeof(double));

ocirc = (fftw_complex*) fftw_malloc(Ny * (Nx / 2 + 1) * sizeof(fftw_complex));

// revv - matrix(c(rev(v), rep(0, C$Kx*(C$Ly-1))), // C$Kx, ncol(C$Cblock));

// revv - rbind(revv, matrix(0, (C$Lx-1), ncol(revv)));

// // mult - fft(C$Cblock * fft(revv), inverse = TRUE);

// // Re((mult/(prod(dim(C$Cblock))))[1:C$Lx,1:C$Ly]);

/* Fill the arrays */ memset(circ, 0, Nx * Ny * sizeof(double));

for (j = 0;

j Ky;

++j) for (i = 0;

i Kx;

++i) circ[i + j*Nx] = v[Kx*Ky - i - j*Kx - 1];

/* Compute the FFT of the reversed vector v */ fftw_execute_dft_r2c(h-r2c_plan, circ, ocirc);

/* Dot-multiply with pre-computed FFT of toeplitz circulant */ for (i = 0;

i Ny * (Nx/2 + 1);

++i) ocirc[i] = ocirc[i] * h-circ_freq[i];

/* Compute the reverse transform to obtain result */ fftw_execute_dft_c2r(h-c2r_plan, ocirc, circ);

/* Cleanup and return */ for (j = 0;

j Ly;

++j) for (i = 0;

i Lx;

++i) out[i + j*Lx] = circ[i + j*Nx] / (Nx*Ny);

fftw_free(circ);

fftw_free(ocirc);

} static void hbhankel_tmatmul(double* out, const double* v, const void* matrix) { const hbhankel_matrix *h = matrix;

R_len_t Nx = h-length.x, Ny = h-length.y;

R_len_t Lx = h-window.x, Ly = h-window.y;

R_len_t Kx = Nx - Lx + 1, Ky = Ny - Ly + 1, i, j;

double *circ;

fftw_complex *ocirc;

/* Allocate needed memory */ circ = (double*) fftw_malloc(Nx * Ny * sizeof(double));

ocirc = (fftw_complex*) fftw_malloc(Ny * (Nx / 2 + 1) * sizeof(fftw_complex));

// revv - matrix(c(rep(0, C$Lx*(C$Ky-1)), rev(v)), // C$Lx, ncol(C$Cblock));

// revv - rbind(matrix(0, (C$Kx-1), ncol(revv)), revv);

// //mult - fft(C$Cblock * fft(revv), inverse = TRUE);

// //Re((mult/(prod(dim(C$Cblock))))[C$Lx:(C$Lx+C$Kx-1), // C$Ly:(C$Ly+C$Ky-1)]);

/* Fill the arrays */ memset(circ, 0, Nx * Ny * sizeof(double));

for (j = 0;

j Ly;

++j) for (i = 0;

i Lx;

++i) circ[(i+Kx-1) + (j+Ky-1)*Nx] = v[Lx*Ly - i - j*Lx - 1];

/* Compute the FFT of the reversed vector v */ fftw_execute_dft_r2c(h-r2c_plan, circ, ocirc);

/* Dot-multiply with pre-computed FFT of toeplitz circulant */ for (i = 0;

i Ny * (Nx/2 + 1);

++i) ocirc[i] = ocirc[i] * h-circ_freq[i];

/* Compute the reverse transform to obtain result */ fftw_execute_dft_c2r(h-c2r_plan, ocirc, circ);

/* Cleanup and return */ for (j = 0;

j Ky;

++j) for (i = 0;

i Kx;

++i) out[i + j*Kx] = circ[(i + Lx - 1) + (j + Ly - 1)*Nx] / (Nx*Ny);

fftw_free(circ);

fftw_free(ocirc);

} static R_INLINE void hbhankelize_fft(double *F, const double *U, const double *V, const hbhankel_matrix* h) { R_len_t Nx = h-length.x, Ny = h-length.y;

R_len_t Lx = h-window.x, Ly = h-window.y;

R_len_t Kx = Nx - Lx + 1, Ky = Ny - Ly + 1;

R_len_t i, j;

R_len_t k, l;

R_len_t wx, dwx, wy, dwy;

double *iU, *iV;

fftw_complex *cU, *cV;

double *wW;

//wW = (double*) fftw_malloc(Nx * Ny * sizeof(double));

/* Allocate needed memory */ iU = (double*) fftw_malloc(Nx * Ny * sizeof(double));

iV = (double*) fftw_malloc(Nx * Ny * sizeof(double));

cU = (fftw_complex*) fftw_malloc(Ny * (Nx / 2 + 1) * sizeof(fftw_complex));

cV = (fftw_complex*) fftw_malloc(Ny * (Nx / 2 + 1) * sizeof(fftw_complex));

/* Fill the arrays */ memset(iU, 0, Nx * Ny * sizeof(double));

for (j = 0;

j Ly;

++j) for (i = 0;

i Lx;

++i) iU[i + j*Nx] = U[i + j*Lx];

memset(iV, 0, Nx * Ny * sizeof(double));

for (j = 0;

j Ky;

++j) for (i = 0;

i Kx;

++i) iV[i + j*Nx] = V[i + j*Kx];

/* Compute the FFTs */ fftw_execute_dft_r2c(h-r2c_plan, iU, cU);

fftw_execute_dft_r2c(h-r2c_plan, iV, cV);

/* Dot-multiply */ for (i = 0;

i Ny * (Nx/2 + 1);

++i) cU[i] = cU[i] * cV[i];

/* Compute the inverse FFT */ fftw_execute_dft_c2r(h-c2r_plan, cU, iU);

/* Form the result */ for (j = 0, wy = 1, dwy = 1;

j Ny;

++j, wy += dwy) { if (j == Ly - 1) dwy--;

if (j == Ky - 1) /* Do not join two ifs! */ dwy--;

for (i = 0, wx = 1, dwx = 1;

i Nx;

++i, wx += dwx) { if (i == Lx - 1) dwx--;

if (i == Kx - 1) dwx--;

F[i+j*Nx] = iU[i+j*Nx] / wx / wy / Nx / Ny;

} } fftw_free(iU);

fftw_free(iV);

fftw_free(cU);

fftw_free(cV);

} static void hbhmat_finalizer(SEXP ptr) { ext_matrix *e;

hbhankel_matrix *h;

if (TYPEOF(ptr) != EXTPTRSXP) return;

e = R_ExternalPtrAddr(ptr);

if (!e) return;

if (strcmp(e-type, "hbhankel matrix")) return;

h = e-matrix;

free_circulant(h);

Free(h);

Free(e);

R_ClearExternalPtr(ptr);

} SEXP initialize_hbhmat(SEXP F, SEXP windowx, SEXP windowy) { R_len_t Nx, Ny, Lx, Ly;

hbhankel_matrix *h;

ext_matrix *e;

SEXP hbhmat;

int *dimF = INTEGER(getAttrib(F, R_DimSymbol));

Nx = dimF[0];

Ny = dimF[1];

Lx = INTEGER(windowx)[0];

Ly = INTEGER(windowy)[0];

/* Allocate memory */ e = Calloc(1, ext_matrix);

e-type = "hbhankel matrix";

e-mulfn = hbhankel_matmul;

e-tmulfn = hbhankel_tmatmul;

e-ncol = hbhankel_ncol;

e-nrow = hbhankel_nrow;

/* Build toeplitz circulants for hankel matrix */ h = Calloc(1, hbhankel_matrix);

initialize_circulant(h, REAL(F), Nx, Ny, Lx, Ly);

e-matrix = h;

/* Make an external pointer envelope */ hbhmat = R_MakeExternalPtr(e, install("external matrix"), R_NilValue);

R_RegisterCFinalizer(hbhmat, hbhmat_finalizer);

return hbhmat;

} SEXP is_hbhmat(SEXP ptr) { SEXP ans = NILSXP, tchk;

ext_matrix *e = NULL;

PROTECT(ans = allocVector(LGLSXP, 1));

LOGICAL(ans)[0] = 1;

/* Object should be external matrix */ PROTECT(tchk = is_extmat(ptr));

/* pointer itself should not be null */ if (LOGICAL(tchk)[0]) { e = R_ExternalPtrAddr(ptr);

if (!e) LOGICAL(ans)[0] = 0;

} else LOGICAL(ans)[0] = 0;

/* finally, type should be ‘hankel matrix’ */ if (LOGICAL(ans)[0] && e && strcmp(e-type, "hbhankel matrix") != 0) LOGICAL(ans)[0] = 0;

UNPROTECT(2);

return ans;

} SEXP hbhankel_rows(SEXP ptr) { SEXP tchk;

SEXP ans = NILSXP;

/* Perform a type checking */ PROTECT(tchk = is_hbhmat(ptr));

if (LOGICAL(tchk)[0]) { ext_matrix *e = R_ExternalPtrAddr(ptr);

PROTECT(ans = allocVector(INTSXP, 1));

INTEGER(ans)[0] = hbhankel_nrow(e-matrix);

UNPROTECT(1);

} else error("pointer provided is not a hankel block-hankel matrix");

UNPROTECT(1);

return ans;

} SEXP hbhankel_cols(SEXP ptr) { SEXP tchk;

SEXP ans = NILSXP;

/* Perform a type checking */ PROTECT(tchk = is_hbhmat(ptr));

if (LOGICAL(tchk)[0]) { ext_matrix *e = R_ExternalPtrAddr(ptr);

PROTECT(ans = allocVector(INTSXP, 1));

INTEGER(ans)[0] = hbhankel_ncol(e-matrix);

UNPROTECT(1);

} else error("pointer provided is not a hankel block hankel matrix");

UNPROTECT(1);

return ans;

} SEXP hbhmatmul(SEXP hmat, SEXP v, SEXP transposed) { SEXP Y = NILSXP, tchk;

/* Perform a type checking */ PROTECT(tchk = is_hbhmat(hmat));

if (LOGICAL(tchk)[0]) { R_len_t K, L;

ext_matrix *e;

hbhankel_matrix *h;

/* Grab needed data */ e = R_ExternalPtrAddr(hmat);

h = e-matrix;

L = (LOGICAL(transposed)[0] ? hbhankel_ncol(h) :

hbhankel_nrow(h));

K = (LOGICAL(transposed)[0] ? hbhankel_nrow(h) :

hbhankel_ncol(h));

/* Check agains absurd values of inputs */ if (K != length(v)) error("invalid length of input vector ’v’");

/* Allocate output buffer */ PROTECT(Y = allocVector(REALSXP, L));

/* Calculate the product */ if (LOGICAL(transposed)[0]) hbhankel_tmatmul(REAL(Y), REAL(v), h);

else hbhankel_matmul(REAL(Y), REAL(v), h);

UNPROTECT(1);

} else error("pointer provided is not a hankel block-hankel matrix");

UNPROTECT(1);

return Y;

} SEXP hbhankelize_one_fft(SEXP U, SEXP V, SEXP hmat) { SEXP F = NILSXP, tchk;

/* Perform a type checking */ PROTECT(tchk = is_hbhmat(hmat));

if (LOGICAL(tchk)[0]) { ext_matrix *e;

hbhankel_matrix *h;

double *rU = REAL(U), *rV = REAL(V), *rF;

/* Grab needed data */ e = R_ExternalPtrAddr(hmat);

h = e-matrix;

/* Allocate buffer for output */ PROTECT(F = allocVector(REALSXP, h-length.x * h-length.y));

rF = REAL(F);

/* Perform the actual hankelization */ hbhankelize_fft(rF, rU, rV, h);

UNPROTECT(1);

} else error("pointer provided is not a hankel block-hankel matrix");

UNPROTECT(1);

return F;

}

Pages:     | 1 |   ...   | 2 | 3 ||
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.