Цифровая обработка многомерных сигналов

Цифровая обработка многомерных сигналов

ТРИ ЭТАЛОННЫХ КРИТЕРИЯ ОЦЕНКИ КАЧЕСТВА СЖАТИЯ ИЗОБРАЖЕНИЙ

Илюшкина Н.С., Чобану М.К.

Московский энергетический институт (ТУ)

Наиболее распространенный объективный критерий оценки СКО (или MSE) является весьма ненадежным, так как он плохо соответствует системе визуального восприятия человека (HVS). Так же следует заметить, что величина MSE в отдельных случаях может незначительно изменяться при существенном ухудшении качества изображения, а, следовательно, она также как и PSNR не может быть взята за основу в построении оптимальных с визуальной точки зрения систем восстановления изображений.

Проведенные нами исследования показали, что система человеческого восприятия (Human visual system HVS) менее восприимчива к искажениям на низких частотах, чем в высокочастотной области. Это было основной причиной учитывать HVS при разработке популярного стандарта сжатия JPEG. В частности, специально спроектированная таблица корректирующих величин, которая применяется для квантования коэффициентов ДКП в блоках 8×8. Согласно этой таблице высокочастотные коэффициенты квантуются с большим шагом квантования, что приводит к лучшему качеству восстановленного изображения, чем если применять для всех коэффициентов одинаковое квантование.

Выражение для MSE с плавающим окном 8×8 с ДКП, принимая в рассмотрение HVS, примет вид:

(1), где Х_ij коэффициенты ДКП блока 8×8, для которого координаты верхнего левого элемента равны i и j, X_ij^e коэффициенты ДКП соответствующего окна оригинального изображения, Т_с матрица корректирующих факторов. В качестве такой матрицы мы используем таблицу квантования стандарта JPEG.

Таблица 1. Матрица Тс корректирующих факторов для вычисления MSE_HVS

1.6084	2.3396	2.5735	1.6084	1.0723	0.6434	0.5046	0.4219
2.1446	2.1446	1.8382	1.3545	0.9898	0.4437	0.4289	0.4679
1.8382	1.9796	1.6084	1.0723	0.6434	0.4515	0.3730	0.4596
1.8382	1.5138	1.1698	0.8874	0.5046	0.2958	0.3217	0.4151
1.4297	1.1698	0.6955	0.4596	0.3785	0.2361	0.2499	0.3342
1.0723	0.7353	0.4679	0.4021	0.3177	0.2475	0.2277	0.2797
0.5252	0.4021	0.3299	0.2958	0.2499	0.2127	0.2145	0.2548
0.3574	0.2797	0.2709	0.2626	0.2298	0.2574	0.2499	0.2600

Матрица Тс обладает следующими двумя свойствами:

1) Соотношения ее коэффициентов обратно пропорциональны соответствующим коэффициентам таблицы квантования JPEG.

2) Суммарный корректирующий фактор матрицы равен единице. Это означает, что в случае однородного распределения искажений по всем пространственным частотам, значения MSE_HVS совпадает с MSE.

Получив формулу для MSE можно записать выражение для PSNR, учитывающего HVS:

Рассмотрим две последовательности {x_i} и {y_i} =1,2,...,N, соответствующие эталонному и восстановленному изображениям, тогда универсальный индекс качества будет определяться комбинацией статистических характеристик (математического ожидания, дисперсии и корреляционной функции) соответствующих сигналов по следующей формуле: (2), где

Значение критерия UQI изменяется в интервале [−1, 1], при этом UQI = 1 соответствует наилучшему качеству изображения. Максимальное значение UQI = 1 получается только тогда, когда x_i= y_i на протяжении всего сигнала (i =1,2,…,N) , а минимальное значение UQI = -1 , в случае если y_i=2-x_i, i=1,2,…,N. Определённый таким образом индекс качества учитывает в себе три искажающих фактора: степень коррелированности отчетов двух сигналов, искажение яркости и искажение контрастности относительно оригинального сигнала. В результате этого исходное выражение может быть представлено в виде произведения трёх множителей:

^(3). Так как в общем случае цифровой сигнал, соответствующий некоторому изображению не является стационарным, необходимо выделить в изображении отдельные локальные области размером A×B, в пределах которых сигнал можно считать стационарным и вычислять статистические характеристики внутри этой области, а уже на основании их вычислять значение UQI_i, соответствующее выделенному блоку. Таким образом, UQI, характеризующий изображение в целом будет определяться как среднее арифметическое значений UQI_i на протяжении всего изображения, то есть , где M - число блоков, внутри которых вычислялся UQI.

Любую вещественную матрицу М можно представить следующим образом: , где U и V единичные матрицы, матрица S размером m×n, диагональные элементы которой не отрицательны и называются сингулярными числами, а вне диагонали равны 0. Такое разложение называется сингулярным разложением матрицы М или разложения по сингулярным числам матрицы.

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

Предложенный графический критерий – это двумерная величина, которая показывает разницу между сингулярными числами исходного и искаженного блока изображения: , где s_i сингулярные числа блока исходного изображения, ŝ_i сингулярные числа блока искаженного изображения, n размер блока. Если размер блока k, получим (k/n)×(k/n) блоков. Множество полученных разностей при графическом отображении представляет "карту искажений".

Для 8-битного черного изображения размером n со значениями пикселей равными 0, сингулярные числа равны 0, для белого изображения размером n со значениями пикселей равными 255, одно сингулярное число равно n×255, а n-1 чисел равны 0.

На рисунке представлены результаты для изображений Monarch и Stream Bridge. Карта искажений, которая отображает степень искажения, тип искажения и ошибку распределения, представляет собой черно-белое изображение. Эта карта получается в результате отображения значений D_i в диапазоне [0,255]. При этом размеры карты искажения уменьшаются на величину блока. Более темным областям соответствует минимальное искажение, а светлым большая степень искажений.

Рис.1. Искаженные изображения и их карты искажений.

Критерий в численном виде вычисляется из графического представления. Рассчитывается общая ошибка, представленная одним числом, в зависимости от типа искажения:

(4), где D_mid среднее значение множества D_i, k размер изображения и n размер блока.

Для сравнения предложенных критериев были проведены субъективные тесты. Искажения были получены для фильтров: 2bior4.4, K5N3, 2bern1.1.2, 2.bern2.2.2, 9db7; алгоритмов: SPIHT и SQP, для различных битрейтов (рис. 2-4).

Рис. 2 Критерий PSNRHVS	Рис. 3 Критерий UQI	Рис. 4 Критерий M-SVD

В тестировании приняло участие 60 человек. Им были представлены пары изображений – оригинал и искаженное изображение, и 10-ти бальная шкала для оценки качества искаженного изображения. 10 соответствовало лучшему качеству, а 1 – худшему. Большинство наблюдателей студенты университета. Размер плавающего окна для каждого метода 8×8. На рис. 2-4 приведены графики зависимости субъективных оценок (MOS) от различных объективных критериев.

В целом в работе критерий PSNR_HVS показал лучшие результаты. Это связано с тем, что в его разработку легло HVS, а именно была использована таблица корректирующих факторов, которая дает представление о восприятии человеком различных частот. Для UQI и M-SVD результаты были получены хуже. Это объясняется тем, что данный критерий не был основан на какой-либо модели HVS. Данные результаты были достигнуты из-за реальной возможности измерения потери качества в процессе искажения изображения при его сжатии. Отличие от PSNR в том, что они чувствительны к энергии ошибок, а не к реальной потере качества.

Литература

1. 
   Egiazarian K., Jaakko A., Ponomarenko N., Lukin V., Battisti F., Carli M., " Two new full-reference quality metrics based on HVS". Proc. 2^nd Intern. Workshop on Video Processing and Quality Metrics for Consumer Electronics, VPQM 2006, Scottsdale, Arizona, USA, 22-24 January 2006, 4p.
   
2. 
   Z. Wang, A. Bovik, “A universal image quality index,” IEEE Signal Proc. Let., 9, 2002, pp. 81–84
   
3. 
   Shnayderman A., Gusev A.,Eskicioglu A., "An SVD-Based Grayscale Image Quality Measure
   

for Local and Global Assessment" IEEE Trans. on Imag. Proc., vol. 15, no. 2, Feb., 2006


THREE FULL-REFERENCE QUALITY MEASURES FOR COMPRESSED IMAGES

Ilushkina N., Tchobanou M.

MPEI (TU)

In this paper three image quality measures are presented. All measures are "full-reference", i.e. they require both the original (reference) image and the distorted one. The traditional Universal Quality Index (UQI) and several recent measures are presented. One of them based on the Peak-Signal-to-Noise Ratio (PSNR) modified to take into account the Human Visual System (HVS). Many studies confirm that the HVS is more sensitive to low frequency distortions rather than to high frequency ones. The other measure is based on singular value decomposition. It can be used as a graphical or as a scalar measure. These measures were used to evaluate wavelet transformation based compression algorithms. Subjective tests were conducted to estimate the performance of the measures. Distorted images were obtained using two algorithms with different filters and bit rates. The analysis of the results shows that PSNR-HVS provides better correlation to Mean Observer Score (MOS) then the other measures.


APPLICATION OF RADIAL BASIS FUNCTION CONCEPT TO MULTIDIMENSIONAL INTERPOLATION PROBLEMS

Tchobanou M.^*, Yamada I.<sup>**

*</sup>Moscow Power Engineering Institute (Technical University), Russia

^**Tokyo Institute of Technology, Japan

In lots of practical applications over a wide field of study one often faces the problem of reconstructing an unknown function f from a finite set of discrete data. These data consist of data sites X = {x₁, . . . , x_N} and data values f_j = f(x_j), 1 ≤ j ≤ N, and the reconstruction has to approximate the data values at the data sites. In other words, a function s is sought that either interpolates the data, i.e. that satisfies s(x_j) = f_j, 1 ≤ j ≤ N, or at least approximates the data, s(x_j) ≈ f_j . The latter case is in particular important if the data contain noise.

Multi-dimensional data, such as images, video signals, usually have high local correlation in all directions. Finding a suitable interpolation algorithm that uses this fact is challenging. The simplest method is to apply a separable algorithm (e.g. by using tensor product functions), which does not address this correlation. However, their implementation tends to be fast, and the mathematics are straightforward. Unfortunately, these decompositions introduce preferred (vertical and horizontal) directions and create a “diagonal” cross term that does not have a straightforward interpretation. This has motivated researchers to design better nonseparable wavelet transforms [1]. One option is to privilege angular selectivity. Numerous directional wavelet transforms, both frames (i.e., redundant) and bases (i.e., nonredundant), have been proposed. Such representations can serve to sparsely represent essential image features such as edges combined with their orientation. Another interesting option, which has received less attention, is to emphasize isotropy. A strong motivation for this kind of design is that many standard image processing algorithms exploit the rotation-invariant properties of filters such as the Gaussian and Laplacian.

Radial basis functions are well-known as traditional and powerful tools for multivariate interpolation from data [2-4] and might be successfully used in different image and video processing applications. Just very recently, radial basis functions have gained enormous popularity in meshfree methods for partial differential equations (PDEs). The theory includes meshfree Galerkin methods, collocation methods, and multilevel schemes. First applications of radial basis functions in computational fluid dynamics are dating back to Kansa. There is nowadays a vast amount of literature on the subject, including using radial basis functions for solving PDEs [5].

The advantage in using radial basis functions is generally that they provide interpolants irrespective of the geometry of the centres X and for any dimension N, whereas for instance polynomial interpolation or even piecewise polynomials are no longer available or feasible in practical applications when N is very large or the geometry of X is complicated.

Suppose a data vector of function values, sampled from an unknown function at a finite point set , d ≥1, is given. Data interpolation requires computing a suitable interpolant satisfying, i.e. for all 1 ≤ j≤ N. Then, the radial basis function interpolation scheme works with a fixed radial function , and the interpolant s is assumed to have the form , (1), where is the Euclidean norm on and denotes the linear space containing all real-valued polynomials in d variables of degree at most m -1, where is said to be the order of the basis function .

Classical choices for radial basis functions , along with their order m, are shown in Table 1, where for any , the symbol denotes as usual the smallest integer greater than or equal to x, and denotes the largest integer less than or equal to x.

Among the most popular radial basis functions are the polyharmonic splines, which are the unique solution of minimizing the functional (2) among all functions in the Sobolev space interpolating on . This class of radial basis functions includes the thin plate splines, where and , which are particularly suited for interpolation from planar scattered data. Further commonly used radial basis functions are given by the Gaussians, , the multiquadrics, of order , and the inverse multiquadrics, , where . Table 1. Radial basis functions

Radial basis function		Parameters	Order
Polyharmonic Splines
Gaussians
Multiquadrics
Inverse Multiquadrics

More recent developments [6-7] have provided compactly supported radial basis functions. In this case, we have for their order, and so the polynomial part in (1) is omitted. While the radial basis functions in Table 1 can be used in arbitrary space dimension d, the selection of one suitable compactly supported depends on d, see Table 2. Since the dimension d is known beforehand, this is no severe restriction. The basics for Wendland’s functions [6] are of the form

, (3), where is a specific univariate polynomial of degree , and so the support of is normalized to the unit interval [0, 1]. Due to Wendland’s construction in [6], the basis function has derivatives up to order 2k. Possible choices for are given in the Table 2, where the symbol denotes equality up to a positive factor, and the truncated power function is given by , for x > 0, and , for x ≤ 0. By their construction, Wendland’s radial basis functions are positive definite on .

DEFINITION. A continuous radial function is said to be positive definite on ,or , iff for any finite set , of pairwise distinct points the matrix is positive definite.

Due to the construction in [7], there exists, for any space dimension d, a positive definite and compactly supported of the form (3). Remarkably enough, Wendland showed that any basis function , constructed in [6] (such as any in Table 2), has minimal degree among all positive definite functions of the form (3). Moreover, by these properties, in (3) is unique up to a positive constant.

The problem of the well-posedness of the interpolation problem (1) is very important. One needs to distinguish the case, where m = 0 from the one where m > 0. First suppose m = 0 for the order of the basis function , such as for the Gaussians, the inverse multiquadrics (in Table 1) and Wendland’s functions (in Table 2). In this case, the interpolant s in (1) has the form

(4). By requiring the n interpolation conditions in (1), the computation of the unknown coefficients of s in (4) amounts to solving the linear equation system (5)

Recall that according to Definition 1, the matrix in (5) is guaranteed to be positive definite, provided that . In this case, the system (5) has a unique solution. This in turn implies the well-posedness of the given interpolation problem already.

Table 2. Compactly supported radial basis functions (Wendland [6])

Dimension d	Radial basis function	Smoothness 2k

Let us turn to the case, where m > 0 for the order of . In this case, the interpolant s in (1) contains a nontrivial polynomial part, yielding q additional degrees of freedom, where is the dimension of the polynomial space . These additional degrees of freedom are usually eliminated by requiring the q moment conditions , for all . (6)

This amounts to solving the linear system (7) where , and for the coefficients of the polynomial part in (1).

One should mention, that classical radial basis functions are poorly conditioned, thus making it very difficult to implement when there are many data points as is typical in signal processing. In addition, the method is computationally very expensive. Non-separable fractional polyharmonic B-splines basis functions are localized versions of radial basis functions, and thus span the same space. If to work out the solution for the fractional case, one may obtain a suitable algorithm for fractal-like signals.

One of possible ways to achieve it might be building of “bell”-shaped splines or B-splines. As it is well known there are 1D B-splines that meet the requirements: a) to get a multiresolution analysis; b) to have an approximation of order L.There are also some results, connected with Laplacean operator approximation (localization filter) in M-D case – namely polyharmonic B-splines [8]. Another ways of achieving the goal might be: a) application of another approximation for M-D Laplacean operator (localization filter); b) building of complex valued polyharmonic B-splines.

The fact that one deals with a uniform grid allows one to develop a fast Fourier-based filtering algorithm. In addition the algorithm can work in any number of dimensions without any special modifications. Thus, one can solve the statistics formulation (Duchon’s splines) with signal-processing techniques (Fourier filtering), achieving a fast algorithm [9].

REFERENCES

[1] TCHOBANOU M., Multi-dimensional multirate systems and multidimensional wavelet functions. Part II. Synthesis. Vestnik MPEI, (3):69–78, 2003.

[2] BUHMANN M.D., Radial basis functions, Acta Numerica (2000), 1–38.

[3] POWELL M.J.D., The theory of radial basis function approximation in 1990, in: “Advances in numerical analysis II: wavelets, subdivision and radial basis functions”, (Ed. Light W.A.), Clarendon Press, Oxford 1992, 105–210.

[4] SCHABACK R., Multivariate interpolation and approximation by translates of a basis function, in: “Approximation theory VIII, Vol. 1: approximation and interpolation”, (Eds. Chui C.K. and Schumaker L.L.), World Scientific, Singapore 1995, 491–514.

[5] Special issue on radial basis functions and partial differential equations, (Eds. Kansa E.J. and Hon Y.C.), Comput. Math. Appl. 43 (2002).

[6] WENDLAND H., Piecewise polynomial, positive definite and compactly sup­ported radial functions of minimal degree, Advances in Comp. Math. 4 (1995), 389–396.

[7] WU Z., Multivariate compactly supported positive definite radial functions, Ad­vances in Comp. Math. 4 (1995), 283–292.

[8] VAN DE VILLE D., BLU T., UNSER M., Isotropic polyharmonic B-splines: scaling functions and wavelets. IEEE Trans. Image Proc., 14(11):1798-1813, Nov 2005.

[9] ABE Y., IIGUNI Y., Fast computation of RBF coefficients using FFT. Signal Processing 86 (2006) 3264–3274



КОЛЕБАНИЯ В АВТОНОМНЫХ ДВУМЕРНЫХ РЕКУРСИВНЫХ ЦИФРОВЫХ СИСТЕМАХ ВТОРОГО ПОРЯДКА С ТРЕМЯ УРОВНЯМИ КВАНТОВАНИЯ

Лебедев М.В.

Ярославский государственный университет им. П.Г. Демидова

Двумерные цифровые системы рекурсивного и нерекурсивного типов широко используются для обработки многомерных цифровых сигналов, статических и динамических изображений [1]. На их основе создаются двумерные цифровые фильтры. Благодаря простоте исполнения и возможности работать в реальном масштабе времени, двумерные цифровые фильтры малых порядков можно использовать в качестве базовых составляющих более сложных устройств цифровой обработки двумерных сигналов и изображений.

Принципиальным отличием цифровых фильтров от аналоговых является обусловленная ограниченным числом используемых двоичных разрядов конечная точность выполнения арифметических операций и задания коэффициентов фильтра, что вызывает специфические ошибки квантования. Исследование процессов в фильтрах может быть осуществлено с помощью линейной статистической модели ошибок, если последовательность ошибок квантования является совокупностью выборок стационарного случайного процесса, если она не коррелированна с последовательностью точных значений сигнала, а сами значения ошибки не коррелированны между собой (представляют собой белый шум) и при этом распределение вероятностей ошибки равномерно во всем диапазоне ошибок квантования [2]. Эти условия нарушаются при малом количестве разрядов (и связанном с ним количестве уровней квантования), а также при воздействии, например, постоянного сигнала или синусоидального сигнала, дискретизированного с частотой, рационально кратной частоте синусоиды. В этих случаях задача исследования процессов в фильтрах очень сложна и является существенно нелинейной [3].

В общем случае движения в автономных двумерных рекурсивных цифровых фильтрах второго порядка с симметричными коэффициентами описывается нелинейным разностным уравнением

(1). Здесь m и n дискретные переменные, принимающие значения от –2 до бесконечности; a, b и c - независимые коэффициенты фильтра. Функция f описывает нелинейные свойства сумматора.

Целью работы является исследование колебаний в автономных двумерных рекурсивных цифровых системах второго порядка с симметричными коэффициентами, использующих представление чисел в прямом коде с тремя уровнями квантования. Следует заметить, что данные колебания являются паразитными для автономных двумерных рекурсивных цифровых фильтров, основанных на исследуемых системах. Полагается, что коэффициенты системы a, b и c задаются без ошибки, квантование осуществляется с округлением, переменные представляются в форме чисел, выровненных справа (т.е. в виде целых чисел), а сумматор без учета квантования имеет характеристику с насыщением.

В качестве характеристики сумматора рассматривается функция (2)

Участок характеристики, соответствующий , назовем зоной I; - зоной II; - зоной III.

Для исследования систем с нелинейным сумматором и квантователями используется детерминированный подход. В частности, он применялся при изучении нелинейных свойств одномерных цифровых фильтров [4-5]. Детерминированный подход используется для решения задач, связанных с изучением условий зарождения предельных циклов разных периодов в результате нелинейной характеристики сумматора в двумерных системах [6].

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

Следует заметить, что любое начальное условие для двумерной системы второго порядка представляет собой две бесконечные последовательности (, и ,), поэтому перебрать все возможные начальные условия при анализе видов движений в такой системе нереально. В связи с этим, при исследовании двумерных систем с произвольными начальными условиями предлагается находить условия возникновения на выходе фильтра только типов движений (каждому из них может соответствовать бесконечное число видов движений), например, двумерных предельных циклов (ДПЦ) [6]. В этом случае достаточно основываться на определении этих циклов и аналитическом виде функции нелинейности.

Рассмотрим это более подробно на примере нахождения областей существования таких типов движений, как ДПЦ с периодами 1х0, 0х1 и 1х1. Первая цифра периода ДПЦ обозначает период по переменной , вторая – по переменной . В случае, если одна из цифр равна нулю, то сигнал может быть как периодическим по соответствующей переменной, так и непериодическим. Предполагается, что отсутствие сигнала на выходе фильтра является частным случаем циклов данных периодов.

Исходя из уравнения (1) и определения данного ДПЦ, найдем область существования цикла периода 1х0

Согласно (2) отсчет может принадлежать одной из зон функции нелинейности. Рассмотрим каждый случай отдельно. Пусть зоне I, т.е.

что равносильно выражению

Следовательно, согласно (2) имеем

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

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

Случаю зоне II соответствует набор условий на коэффициенты фильтра:

На плоскости коэффициентов фильтра (,) при фиксированном коэффициенте данным ограничениям соответствует заштрихованная область 1 на рис. 1.

Вследствие симметричности функции (2) случаю зоне III соответствует такой же набор неравенств, как и в случае, зоне I. В связи с этим разбиение пространства коэффициентов фильтра сохраняет вид, показанный на рис. 1.

Вследствие симметричности коэффициентов и в уравнении (1) относительно переменных и , замена соответствует замене . Таким образом, для ДПЦ с периодом 0х1 остается справедливым полученный набор неравенств, с учетом соответствующих переобозначений. На плоскости коэффициентов (,) при данным циклам на выходе фильтра соответствует заштрихованная область 2 на рис. 1.

Для существования невырожденных циклов 1х0 и 0х1 необходимо, чтобы минимум один из отчетов последовательностей , , и был ненулевым (в противном случае, сигнал на выходе фильтра будет отсутствовать).

Найдем условия на коэффициенты фильтра, соответствующие ДПЦ периода 1х1. Согласно определению данного ДПЦ и уравнению (1), справедлива система уравнений

Отсчет может принадлежать одной из трех зон характеристики сумматора. Рассмотрим каждый случай отдельно. Пусть зоне I, тогда _, что равносильно набору условий на коэффициенты фильтра.

Случаю зоне II соответствует набор условий на коэффициенты фильтра:

На плоскости (,) им соответствует заштрихованная область 3, показанная на рис. 1.

Вследствие симметричности функции нелинейности (2), набор условий и разбиение плоскости (,) для случая зоне III сохраняются неизменными.

При коэффициенте разбиение плоскости (, ) также имеет вид, показанный на рис. 1. Таким образом, представленная на рис. 1 область 3 соответствует циклам периода 1х1 для , принадлежащего любой из зон функции нелинейности. В отличие от циклов других периодов для ДПЦ с таким периодом любой ненулевой отсчет последовательностей начальных условий приводит к образованию диагонали аналогичных отсчетов на выходе фильтра.

Следует заметить, что если один из коэффициентов фильтра больше единицы, а два других лежат в интервале между –0.25 и 0.25, на выходе возникает ДПЦ с периодом, зависящим от того, какой именно из коэффициентов, превосходит единицу. В случае это ДПЦ периода 1х0, если , то имеем цикл периода 0х1, условию соответствует ДПЦ периода 1х1.

Рис 1. Бифуркационная диаграмма двумерного рекурсивного цифрового фильтра второго порядка с симметричными коэффициентами

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

Литература

1. 
   Даджион Д., Мерсеро Р. Цифровая обработка многомерных сигналов. - М.: Мир, 1988. - 488 с.
   
2. 
   Каппелини В., Константинидис А. Дж., Эмилиани П. Цифровые фильтры и их применение. - М.: Энергоатомиздат, 1983. - 360 с/
   
3. 
   T. Bose, D. Brown, “Zero-input limit cycles due to rounding in digital filters”, IEEE Trans. Circuits & Syst., V. 6, 1989. P. 931-933.
   
4. 
   Брюханов Ю.А. Эффекты квантования в цифровых рекурсивных фильтрах первого порядка с усечением по величине // Известия вузов. Прикладная нелинейная динамика. Т. 10, № 6. 2002. С. 35-41.
   
5. 
   Брюханов Ю.А. Колебания в цифровых рекурсивных фильтрах первого порядка с представлением чисел в дополнительном коде с округлением // Известия вузов. Прикладная нелинейная динамика. Т. 12, № 1-2. 2004. С. 10-17.
   
6. 
   D.V. Rudyh, M.V. Lebedev, V.V. Kryashchev, and A.L. Priorov. Investigation of the two-dimensional first-order recursive digital filters with saturation nonlinearity // Proc. of the 11-th Workshop on “Nonlinear Dynamics of Eletronic Systems” (NDES'2003), Switzerland, 2003. P. 213-216.
   


OSCILLATIONS IN AUTONOMUS TWO-DIMENSIONAL RECURSIVE DIGITAL SYSTEMS OF SECOND ORDER WITH THREE LEVELS QUANTISATION

Lebedev M.

Yaroslavl state university by P.G. Demidov
Two-dimensional digital systems of recursive and non-recursive types are widely used for processing multivariate digital signals, static and dynamic images. On their basis two-dimensional digital filters are created. Owing to simplicity of execution and an opportunity to work in real time, two-dimensional digital filters of small orders can be used as base making more complex devices of digital processing two-dimensional signals and images.

Generally movements in independent two-dimensional recursive digital filters of the second order with symmetric coefficients it is described nonlinear by the equation

The purpose of work is research of fluctuations in independent two-dimensional recursive digital systems of the second order with the symmetric coefficients, numbers using representation in a direct code with three levels of quantization. Dynamics is considered on an example of a finding of areas of existence of such types of movements with the periods 1x0, 0x1 and 1x1. On pic.1 the allocated areas correspond with the periods 1x0, 0x1 and 1x1.

Fig. 1. Bifurcation diagram of the two-dimensional recursive digital system of the second order with symmetrical coefficients

Thus, the technique is presented, allowing finding areas in space of coefficients of the filter, corresponding the set types of movements. With its help areas of two-dimensional limiting cycles of the various periods are certain. It is necessary to notice, that the given types of movements are parasitic for the two-dimensional recursive digital filters based on investigated systems. Analytical conditions on coefficients of the two-dimensional recursive digital filter of the second order with symmetric are found. The diagram of system is constructed. Theoretical results are confirmed by computer modeling.



Алгоритм оценки параметров и экстраполяции двухмерного сигнала имеющего гармоническую составляющую

Карамов С.В.

ЗАО «НТЦ ЭЛИНС»,

Россия, 124460, Зеленоград, Панфиловский проспект, д.4, стр.1, тел.(495)5320151, e-mail: sergkar@elins.ru

Рассматривается задача оценки параметров и экстраполяции двумерного дискретного сигнала имеющего гармоническую составляющую. На плоскости подобный сигнал отображается в виде дискретной траектории, состоящей из поступательной и вращательной составляющих. По совокупности входных отсчетов представляющих двумерные координаты, необходимо оценить параметры сигнала и осуществить экстраполяцию траектории при возможном отсутствии части измеренных координат. Данная задача появляется, к примеру, при слежении и управлении перемещающимися в пространстве вращающимися летательными аппаратами. Физически измерение координат возможно с помощью оптико-электронных систем, радиолокации, либо иными методами. В параметрическом виде модель измеренного двумерного сигнала представляется в виде уравнений: (1)? где: – координаты поступательной составляющей (центра вращения); – радиус вращения; – фаза вращения; – шумы измерения. Нестационарные параметры сигнала (1) в общем случае могут изменяться произвольно. Существуют несколько подходов к оценке параметров и прогнозированию подобного сигнала [1-3].

Предлагается новый метод, основанный на нелинейной локальной аппроксимации двумерной траектории некоторой моделью. Для определения параметров модели в текущем временном отсчете, будем аппроксимировать траекторию на локальном участке длительностью T, содержащем предыдущих отсчетов. Далее рассчитанные параметры поступают на фильтр Калмана, работающий в режиме прогноза, с выхода которого получаем сглаженные значения параметров. Затем происходит восстановление траектории по модели (1) с подавлением шума и помех.

Используются равноотстоящие отсчеты с интервалом дискретизации . Частота дискретизации сигнала . Введем признак наличия достоверных координат e_k. Отсчеты, содержащие информацию с достоверными координатами, имеют , в противном случае , где – номер отсчета в текущем анализируемом временном интервале. Будем аппроксимировать параметры модели сигнала под измеренный входной сигнал, в предположении квадратичной изменения поступательного движения, линейной аппроксимации изменения радиуса и фазы вращения: . (2)

Вводя обозначения: , , , , , , , , перепишем (2) в виде: . (3)

Для нахождения одиннадцати неизвестных параметров: x_o, v_x, a_x, y_o, v_y, a_y, R₁, ρ₁, R₂, ρ₂, ω, воспользуемся нелинейным методом наименьших квадратов.

Определим расстояние от измеренных координат (x_k,y_k) до точки лежащей на аппроксимирующей кривой для моментов времени как: . Минимизируем по отсчетам функцию: , для чего решим систему из одиннадцати уравнений: . (4)

Из-за громоздкости получаемых выражений, опустим вывод развернутых формул. Поскольку систему (4), в которой искомые параметры связаны нелинейными зависимостями, затруднительно решить аналитически, решим ее итерационно численным методом. Зафиксируем параметр ω, и из первых десяти уравнений получим систему линейных уравнений относительно параметров x_o, v_x, a_x, y_o, v_y, a_y, R₁, ρ₁, R₂, ρ₂. Запишем сокращенную систему в матричном виде: , (5)? где: ,

,, , , , , , , , , , , ,, , , , , , , , , , , , , , , , , , , .

Исключенное из (4) одиннадцатое уравнение представляется:

(6)

Если для : , то вычисление конечных сумм значительно упрощается, используя их представление в явном виде по известным формулам. Вектор параметров из (5) легко находится обычным способом как решение системы линейных уравнений. Среднеквадратическая ошибка аппроксимации χ, рассчитывается по формуле:

Таким образом, задача сводится к минимизации одномерной функции , определенной на отрезке возможных значений параметра ω. Необходимо учитывать, что χ(ω) может иметь несколько локальных минимумов, поэтому необходима процедура поиска глобального минимума на отрезке допустимых частот вращения. Для найденной ω, соответствующей глобальному минимуму χ(ω), рассчитываются оставшиеся параметры x_o, v_x, a_x, y_o, v_y, a_y, R₁, ρ₁, R₂, ρ₂. Откуда выражаются , .

Функция χ(ω) имеет пульсирующий характер с множеством локальных минимумов. Локальные минимумы расположены по обе стороны относительно глобального минимума, достигаемого на частоте ω, и имеют амплитуды, возрастающие по мере приближения к концам диапазона . Аддитивный шум искажает амплитуды локальных минимумов и может уменьшить их число. Период пульсаций локальных минимумов равен (Гц). На диапазоне количество минимумов не превышает значения (n-1).

В случае априорно не известной частоты вращения и направления вращения, для расчета используется весь диапазон. Причем, при использовании в процессе измерений правой декартовой системы координат, соответствует вращению против часовой стрелки, а соответствуют вращению по часовой стрелке. Минимально необходимое количество значений ω при первичном поиске лепестка глобального минимума χ(ω) на отрезке , должно быть не менее 2n. При априорно известном диапазоне частот , минимально необходимое количество значений ω не менее , где int( ) – округление до большего целого. Далее, после первичного поиска, рассчитав ω соответствующую минимуму χ(ω), любым из известных итерационных методов достигаем глобального минимума χ(ω). Необходимо учитывать, что при любых значениях параметров сигнала и величинах шума, χ(ω) остается гладкой функцией, что облегчает нахождение ее экстремального значения. Дополнительно, для ускорения поиска глобального минимума можно использовать невязку μ(ω), получаемую из (6), пересекающую ноль с положительным наклоном в точках минимума χ(ω) а в окрестности глобального минимума имеющую максимальные на всем диапазоне биполярные выбросы.

Рис 1. Ошибка χ(ω) Рис. 2 Невязка μ(ω)

При необходимости ограничиться меньшим числом параметров в модели (2), что часто случается на практике, из H требуется исключить соответствующие строки и столбцы. Использование меньшего числа параметров при аппроксимации повышает устойчивость алгоритма при большом количестве недостоверных отсчетов и значительно понижает число обусловленности матрицы H. При использовании полного набора параметров, необходимо соответственно, больше итераций при поиске ω, поскольку лепесток глобального минимума имеет большую крутизну и характерные высокоамплитудные боковые составляющие.

На рис.1 представлен характер зависимости ошибки χ(ω) при использовании параметра ρ, при f=20 (Гц), f_d=200 (Гц). На рис.2 представлен характер зависимости невязки μ(ω).

Рассчитанные значения локальных параметров имеют методические погрешности вычисления и подвержены во времени низкочастотным эволюциям. В связи с этим, для сглаживания мгновенных значений полученных параметров и их прогнозирования в случае пропадания координат используется фильтр Калмана [1], работающий в режиме прогноза, порядок которого выбирается в соответствии с предполагаемой динамикой сигналов. По полученным на выходе фильтра Калмана сглаженным оценкам параметров , происходит восстановление исходной сглаженной траектории в соответствии с предполагаемой моделью сигналов (1). При пропадании входных координат, восстановление траектории источника происходит по (1) с учетом:, где: – момент пропадания сигнала; – текущее рассчитанное значение фазы вращения в момент пропадания координат.

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

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

Литература

1. 
   Понятский В.М., Карамов С.В. Интеллектуальная система формирования строба в видеокадре при слежении за вращающимся источником полезного излучения // Материалы межрегиональной научно-технической конференции «Интеллектуальные и информационные системы». Тула, 2004. –С 53-56.
   
2. 
   Карамов С.В. Оценка параметров и прогноз движения вращающегося объекта имеющего трохоидальную траекторию по видеоизображению // Труды XVI Международная конференция по компьютерной графике и ее приложениям «Графикон’2006». Новосибирск, 2006.–С 347-350.
   
3. 
   Карамов С.В., Понятский В.М. Двухмерный квазиоптимальный алгоритм оценки параметров траектории движения вращающегося летательного аппарата // Изв. ТулГУ. Сер. «Проблемы проектирования и производства систем и комплексов». Вып. 9 (ч.1). –Тула: ТулГУ, 2006.–С 99-103.
   



ALGORITHM OF THE ESTIMATION OF PARAMETERS AND EXTRAPOLATIONS OF THE two-dimensional SIGNAL HAVING THE harmonic COMPONENT

Karamov S.

Joint-Stock Company «NTC ELINS»,

Stroenie 1, Panfilovskiy prospect 4, 124460, Zelenograd, Russia, tel.(495)5320151, e-mail: [email protected]

The problem of an estimation of parameters and extrapolations of a two-dimensional discrete signal having a harmonic component is considered. On a plane the similar signal is displayed in the form of the discrete trajectory consisting of forward and rotary components. On set of the input samples representing two-dimensional coordinates, it is necessary to estimate parameters of a signal and to carry out extrapolation of a trajectory at possible absence of a part of the measured coordinates. The given problem appears, for example, at tracking and control moving in space rotating flying. Physically measurement of coordinates probably by means of opto-electronic systems, a radar, or other methods.

In a parametrical kind the model of the measured two-dimensional signal is represented in the form of the equations:

(1), where: – coordinates of a forward component (the center of rotation); – radius of rotation; – a phase of rotation; – noise of measurement. Non-stationary parameters of a signal (1) generally can be any. There are some approaches to an estimation of parameters and forecasting of a similar signal. The new method based on local nonlinear approximation of a two-dimensional trajectory is offered. For definition of parameters of model in current samples, the input signal is approximated on the local site containing of the previous readout. Further the calculated parameters act on filter Kalman working in a mode of the forecast from which output it is received the smoothed values of parameters. Then there is a restoration of a trajectory according to model (1) to suppression of noise and handicaps. Equidistant time samples are used. Parabolic approximation of progress, linear approximations of radius and a phase of rotation in model of a kind is supposed: .

For minimization mean-square error of approximation, the iterative numerical method with a finding of a global minimum is offered. Character of function of a error of approximation depending on parameters is analyzed. The offered algorithm allows to carry out with high accuracy an estimation of parameters and extrapolation a signal under condition of incomplete entrance data.

