Теория и методы цифровой обработки сигналов

plot(fe,10*log10(abs(EZ2(1:128))),['K','-'])

hold on

plot(f2,10*log10(abs(Z2(1:128))),['K','--'])

%ПЕРЕКРЕСТНОЕ ИСПОЛЬЗОВАНИЕ КОЭФФИЦИЕНТОВ ОТРАЖЕНИЯ

for i = 2:p+1,

ep = ef(i:N); ef-ошибка предсказания вперед первого канала

em = eb(i-1:N-1);

K(i-1) = 2 * ep' * em / (ep'*ep + em'*em);

a = [a;0] - K(i-1) * [0;flipud(a)];

% ПЕРЕКРЕСТНЫЙ ВВОД КОЭФФИЦИЕНТА ОТРАЖЕНИЯ ВТОРОГО КАНАЛА К2

% ПРИ ФОРМИРОВАНИИ ОШИБКИ ПРЕДСКАЗАНИЯ ef ПЕРВОГО КАНАЛА

for j = N:-1:i,

ef_old = ef(j);

ef(j) = ef(j) - K2(i-1) * eb(j-1);

eb(j) = eb(j-1) - K2(i-1) * ef_old;

end

E(i) = (1 - K(i-1)'*K(i-1)) * E(i-1);

end

figure(4) % Спектральное представление сигналов на выходе решетчатого фильтра с %перекрестными связями

EZ = fft(ef',256);% для 2канала

fe = 256*(0:127)/256;

plot(fe,10*log10(abs(EZ(1:128))),['K','-'])

hold on

plot(f2,10*log10(abs(Z1(1:128))),['K','--'])

Из графика на Рис.4 видно как двухмодовая помеха режектируется адаптивным решетчатым фильтром, а полезный сигнал нет. Аналогичным образом можно показать работу предложенного способа при использовании вобуляции периода повторения импульсов от пачки к пачке. Главный вывод - проведенное исследование в системе MATLAB полностью подтверждает положительный эффект от применения предложенного способа фильтрации дискретных помех по неклассифицированной выборке наблюдений.

Литература

1. 
   Haykin S., Currie B., Kesler S. “Maximum Entropy spectral analysis of radar clutter”, Proc. IEEE Trans.vol.70, September 1982, pp.953-962.
   
2. 
   Barbarossa B., Picardi G. “Predictive Adaptive Moving Target Indicator”, Signal Processing, №10, 1986, pp.83-97.
   
3. 
   Потемкин В.Г. “Справочник по MATLAB” Анализ и обработка данных.
   

http://matlab.exponenta.ru/ml/book2/chapter8/


Adaptive lattice filter for suppression of the discrete correlated interference

Bartenev V.

The construction of the adaptive lattice filter is offered allowing to improve efficiency of suppression of discrete correlated interference. Indicated technical result is achieved due to forming the coefficients of reflection using two different carriers, or two repetition frequencies. The formed coefficients of reflection proposed to use for calculation of errors of prediction to and fro in a such way. Thus coefficients of reflection formed on one carrier or repetition frequency use during filtration on the other carrier or repetition frequency and vice versa. The choice of difference for carriers and repetition frequencies is such that allow for low speed discrete correlated interference effectively suppressed, but for useful signal of high speed target not compensated.

The rightness of the offered adaptive lattice filter is confirmed by the results of its simulation in Matlab.



Одновременное использование нескольких критериев в адаптивных антенных решетках

Джиган В.И.

ГУП НПЦ «ЭЛВИС»

а/я 19, Центральный пр-т, Зеленоград, Москва К-460, 124460. Тел.: +7-095-531-1961. [email protected]. http://multicore.ru

Управление адаптивными антенными решетками (ААР) является одним из применений алгоритмов адаптивной фильтрации. В процессе выполнения этих алгоритмов в диаграмме направленности (ДН) ААР формируются провалы в направлениях источников помех, в результате чего помехи подавляются. Особенностью функционирования ААР в системах связи является частое отсутствие источника требуемого сигнала, которое заменяется знанием направления на источник полезного сигнала или знанием его формы. В первом случае адаптивные алгоритмы строятся на основе линейно-ограниченных критериев (Linearly-Constrained, LC) [1], а во втором – на основе критерия постоянства модуля огибающей сигналов (Constant Modulus, CM) [2]. Однако известно, что LC-алгоритмы чувствительны к когерентным помехам [3], а СМ-алгоритмы ­в ряде случаев могут «цепляться» за такие помехи, ориентируя основной лепесток ДН ААР в направлении источника помехи и подавляя полезный сигнал [4].

В [5] было показано, что если направление на источник полезного сигнала известно, то введение линейных ограничений в адаптивный алгоритм на основе CM-критерия позволяет ААР эффективно функционировать при наличии когерентных помех. Ограничения удерживают основной лепесток ДН ААР в направлении на полезный сигнал в течение всего времени адаптивной фильтрации. Алгоритм [5] c вычислительной сложностью базируется на стратегии градиентного поиска. Известно, что такие алгоритмы обладают более медленной сходимостью и большими значениями остаточных ошибок в установившемся режиме по сравнению с рекурсивными алгоритмами на основе критерия наименьших квадратов (Recursive Least Squares, RLS) с вычислительной сложностью [6]. Кроме того, минимизируемый функционал в алгоритмах на основе СМ-критерия не является квадратичным. Поэтому градиентные алгоритмы часто сходятся к локальным решениям задачи адаптивной фильтрации.

В настоящей работе рассматривается применение LC RLS-алгоритмов на основе обратного QR-разложения в качестве алгоритмов управления СМ ААР.

Адаптивные алгоритмы на основе СМ-критерия применяются для обработки информационных сигналов, характеризуемых постоянным значением модуля огибающей, например, сигналов с Quadrature Phase Shift Keying (QPSK) модуляцией. Каждый информационный символ при такой модуляции обладает свойством . Значение является известным на приемной стороне. СМ-критерий формулируется как , (1), а адаптивные CM-алгоритмы обозначаются как СМ(). Здесь – выходной сигнал антенной решетки (рис. 1), – вектор весовых коэффициентов, – вектор пространственно-временных отсчетов сигналов, – индекс дискретного времени. В [7] было показано, что при и уравнение (1) может быть преобразовано к виду: , (2), где и .

Т


Рис. 1. Адаптивная антенная решетка без требуемого сигнала
аким образом, функционал (2) является квадратичным. Данное допущение позволяет использовать любой RLS-алгоритм с вычислительной сложностью в качестве алгоритма СМ(2,2). В [7] для этой цели использовался RLS-алгоритм на основе леммы об обращении матриц. Однако исследования показали, что такой алгоритм в СМ ААР может быть неустойчивым.

В настоящей работе рассматривается использование в качестве алгоритма управления СМ ААР многоканального LC RLS-алгоритма на основе устойчивой процедуры обратного QR-разложения с преобразованием Хаусхолдера,. Алгоритм (табл. 1), является частным случаем многоканальных алгоритмов [6]. Он минимизирует функционал (2) при условии , где – вектор фазирования ААР в направлении на источник полезного сигнала , а – значение ДН ААР в направлении .

Таблица 1. LC СМ() RLS-алгоритм

	Вычисления		Ссылки
			(1.0)
			(1.1)
			(1.2)
			(1.3)
			(1.4)
			(1.5)
			(1.6)
			(1.7)
			(1.8)
			(1.9)
			(1.10)
			(1.11)
			(1.12)
			(1.13)
			(1.14)
			(1.15)
а)		б)
в)		г)

Рис. 2. Результаты моделирования

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

Для демонстрации эффективности алгоритма (табл. 1) проведено компьютерное моделирование подавления помех (рис. 2) при , , . Сравнивались 3 алгоритма. В качестве полезного сигнала использовался сигнал с модуляцией QPSK-4 при . Уровень полезного сигнала был принят за 0 дБ. Направление его источник было выбрано как , а направления на источники помех – как и . Коррелированная помеха с уровнем  дБ представляла собой копию полезного сигнала, задержанного на половину длительности информационного символа. Уровень некоррелированной помехи (белого шума) равнялся 20 дБ.

Сравниваемый LC RLS-алгоритм (рис. 2а) отличался от алгоритма (табл. 1) тем, что переменная заменялась на и . СМ(2,2) RLS-алгоритм без ограничений инициализировался как (рис. 2б) и как (рис. 2в), где – единичный вектор. Начальное значение вектора весовых коэффициентов ААР на основе LC RLS-алгоритма (рис. 2а и рис. 2г) определялось как . На рис. 2 серым цветом обозначены ДН синфазной антенной решетки, а черным цветом ­– ДН ААР в установившемся режиме.

Из рис. 2а следует, что LC RLS-алгоритм обеспечивает заданное ограничение и подавление некоррелированной помехи. Подавление коррелированной помехи незначительное и ДН ААР «разваливается». Созвездие информационных символов на выходе ААР из-за присутствия неподавленной коррелированной помехи также разваливается (рис. 3а).

СМ(2,2) RLS-алгоритм (рис. 2б) подавляет обе помехи. В процессе адаптации формируется ДН с направлением , т.е. с направлением на максимальный СМ-сигнал. Однако созвездие приобретает постоянный фазовый сдвиг (рис. 3б) .

В силу иной инициализации СМ(2,2) RLS-алгоритм (рис. 2в) в процессе адаптации «цепляется» за коррелированную помеху, в результате чего ДН переориентируется в направлении этой помехи. Некоррелированная помеха и полезный сигнал подавляются. Созвездие также приобретает фазовый сдвиг (рис.3в).

Для компенсации обычно используется фазовая автоподстройка [2]. Эту компенсацию можно также осуществить, задав ограничение . Моделирование подтверждает, что такое ограничения обеспечивает на каждой итерации не только требуемые ориентацию и уровень основного лепестка ДН ААР (рис. 2г), но и правильную ориентацию созвездия информационных символов (рис. 3г), совпадающую с ориентацией алфавита сообщения.

а)	б)	в)	г)

Рис. 3. Созвездия

Таким образом, введение линейных ограничений в СМ(2,2) алгоритм позволяет эффективно бороться с коррелированными помехами, предотвращая явление «захвата» этих помех, и просто компенсирует постоянный фазовый сдвиг в выходном сигнале ААР.

Литература

1. Frost O.L. An algorithm for linearly constrained adaptive array processing // Proceedings of the IEEE. – 1972. – Vol. 60. – №8. – P. 926–935.

2. Treichler J., Larimore M. New processing techniques based on the constant modulus adaptive algorithm // IEEE Trans. Acoustics, Speech and Signal Processing. – 1985. – Vol. 33. – №2. – P. 420–431.

3. Shan T.-J., Kailath T. Adaptive beamforming for coherent signals and interference // IEEE Trans. Acoustics, Speech, and Signal Processing. – 1985. – Vol. 33. – №3. – P. 527–5361.

4. Treichler J., Larimore M. The tone capture properties of CMA-based interference suppressors // IEEE Trans. Acoustics, Speech, and Signal Processing. – 1985. – Vol. 33. – №4. – P. 946–958.

5. Rude M.J. Griffiths L.J. Incorporation of linear constraints into the constant modulus algorithm // Intern. Conf. on Acoustics, Speech and Signal Processing. – 1989. – Vol. 2. – P. 968–971.

6. Джиган В.И. Многоканальные RLS- и быстрые RLS-алгоритмы адаптивной фильтрации // Успехи современной радиоэлектроники. – 2004. – №11. – С. 48–77.

7. Chen Y., Le-Ngoc T., Champagne B., Xu C. Recursive least squares constant modulus algorithm for blind adaptive array // IEEE Trans. Signal Processing. – 2004. – Vol. 52. – №5. – P. 1452–1456.


Adaptive arrays, simultaneously based on a few criteria

Djigan V.

Electronic VLSI Engineering & Embedded Systems (ELVEES) R&D; Center of Microelectronics

POB 19, Centralny Prospect, Zelenograd, Moscow K-460, Russia 124460

Tel.: +7-095-531-1961. E-mail: [email protected] . URL: http://multicore.ru

T

Fig. 1. Simulation results
he paper considers the application of the linearly constrained (LC) Recursive Least Squares (RLS) algorithm based on inverse QR-decomposition with Householder transform in adaptive antenna arrays (AA) for communication systems with Constant Modulus (CM) signals, QPSK modulated, for example. Following [1, 2], it was demonstrated, that the algorithm provides long-time stability and prevents well-known “tone capture” phenomenon in CM adaptive filtering algorithms. The considered algorithm is a modification of the multichannel RLS algorithms, presented in [3].

Based on the algorithm, eight antennas AA efficiently suppresses both the coherent and white-noise interferences at and respectively, see Fig. 1. Besides, signal constellation at AA output has a correct position which is the same as that of information alphabet, Fig. 2d. At the same time, LC RLS algorithm does not suppress the coherent interference. It destroys directional pattern and AA constellation, Fig. 2a. CM RLS algorithm adds a constant phase shift in AA output signal, Fig. 2b and Fig. 2c. In Fig. 2c case, interference is captured and desired signal is cancelled.

a)	b)	c)	d)

Fig. 2.QPSK-4 constellations at AA output

Thus, if desired signal direction is known, the using of the proposed adaptive filtering algorithm provides efficient performance of AA in coherent interference environment and rotates AA output constellation to correct position.

References

1. Rude M.J. Griffiths L.J. Incorporation of linear constraints into the constant modulus algorithm // Intern. Conf. on Acoustics, Speech and Signal Processing. – 1989. – Vol. 2. – P. 968–971.

2. Chen Y., Le-Ngoc T., Champagne B., Xu C. Recursive least squares constant modulus algorithm for blind adaptive array // IEEE Trans. Signal Processing. – 2004. – Vol. 52. – №5. – P. 1452–1456.

3. Djigan V.I. Multichannel RLS and fast RLS adaptive filtering algorithms // Successes of Modern Radioelectronics. – 2004. – №11. – P. 48–77.


ЦИФРОВОЙ СПЕКТРАЛЬНЫЙ АНАЛИЗ ПОЛИГАРМОНИЧЕСКОГО СИГНАЛА

Баранов И.В. Езерский В.В.

ООО “Предприятие Контакт-1”, г. Рязань

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

Целью данной работы является снижение вычислительной сложности метода и оценка достижимых при этом точностных характеристик.

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

Для решения этой задачи часто вычисляется квадрат метрического расстояния между измеренным сигналом и сигналом с предполагаемыми (перебираемыми) параметрами (1)

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

Этот метод требует выполнения большого объёма вычислений. Развитие метода, предложенное в [1], сводится к выполнению трёх этапов. На первом этапе производится грубая оценка искомых значений координат и амплитуд спектральных линий. В качестве первого приближения искомого спектра берутся амплитуды и координаты хорошо разрешенных сигнальных пиков, полученных при Фурье-преобразовании принятого сигнала:

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

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

Шаги (а) и (б) повторяются до тех пор, пока амплитуды выделяемых пиков не станут меньше выбранного порога.

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

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

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

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

Рассмотрим возможность применения процедуры одномерной минимизации вместо трёхмерной.

Представим искомую функцию в виде (3).

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

;	(4)
.	(5)

Теперь процедуру многомерной численной минимизации можно свести к процедуре одномерной минимизации по частоте . После нахождения частоты , соответствующей минимуму невязки, вычисляем искомые амплитуды и фазы синусоид: , (6), (7)

Проверка возможности использования предложенной модернизации метода оценки параметров сигнала проводилась с помощью компьютерного моделирования. Моделирование проводилось в среде MATLAB. Формировалось отсчётов сигнала. Рассматривались два варианта. В первом варианте моделирования производилось оценка частоты одиночного сигнала, число периодов которого в интервале наблюдения изменялось от 0.5 до 4. При этом сильное влияние на погрешность оценки частоты оказывает спектр в отрицательной области частот. Далее смоделированный сигнал обрабатывался в соответствии с алгоритмом, предложенным выше. Моделировался процесс измерения без шумов для определения предельных возможностей метода. Результаты приведены на рис. 1, где по осям отложены частоты и погрешности в относительных единицах (бинах). Погрешность измерения на 5 порядков меньше, чем при использовании метода Фурье.

Второй вариант моделирования включал моделирование двух сигналов. Первый сигнал полезный, частота которого изменялась. Второй сигнал с постоянной частотой и амплитудой на 3 дБ меньше. Частота полезного сигнала изменялась симметрично относительно частоты мешающего в пределах 3-х бин. Результаты измерения частоты полезного сигнала приведены на рис. 2, где показано половина графика, симметричного относительно 20 бин. На рисунке точечной линией показана погрешность измерения методом Фурье с использованием весового окна Блэкмана и сплошной линией показана погрешность предложенного метода. Видно, что для предложенного метода зона повышенной погрешности измерения равна бин в области мешающего сигнала и затем резко падает. У метода Фурье этот диапазон составляет бина.

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

Литература

1. 
   Свентковский Р. А. Сверхразрешение сигналов: возможности, ограничения, неавторегрессионный подход // Радиотехника и электроника. 1998. №3. С. 288-292.
   


DIGITAL SPECTRAL ANALYSIS POLIHARMONICS SIGNAL
Baranov I. Ezerski V.

LTD “Enterprise Kontakt-1”

The Parameters spectral forming poliharmonic of the signal (the amplitudes , frequencies and phases ) are valued on measurements , got with mistake. It Is offered known three stage procedure [1] reduce to two stage, compatible two last stages and executing searching for spectral forming with simultaneous optimization their parameter on base of the average square of the mistake. The second stage, consisting of two steps, is repeated until difference newly found parameter from parameter, got on previous stage, will not become less certain beforehand given values.

Formulas are received for amplitudes quadrature component and each harmonicas, minimizing average square of the mistake under fixed importance of the frequency

;	(1)
.	(2)

These formulas allow on two step to use the procedure to univariate numerical minimization on frequency instead of three-dimensional minimization. After finding of the frequency , corresponding to minimum discrepancy, are calculated sought amplitude and phase sine wave.

By means of computer modeling is organized checking the possibilities of the offered method and is proved his efficiency

The literature

1. Sventkovskiy R.A. Sverhrazreshenie of signals: opportunities, restrictions, not autoregression the approach // the Radio engineering and electronics. 1998.№3. With. 288-292.





Алгоритм цифровой фильтрации сигналов на основе методов сплайн аппроксимации

Просочкин А.С.

Филиал «Восход» МАИ

468320, Байконур, пр. Гагарина, 5, [email protected]

Для приближения, восстановления и фильтрации сигналов, заданных в виде функций времени f(t) на равномерной сетке t_i (i=−∞, …, −2, −1, 0, 1, 2, …,+∞, Δt=t_i – t_i+1 = const) значениями f(t_i) с некоторой погрешностью, широко используются методы полиномиальной сплайн-аппроксимации [1]. Распространённым способом построения полиномиальных сплайнов произвольной степени m является их представление через базисные B-сплайны в виде [1] , (1), где C_i - коэффициенты аппроксимации; - базисные B-сплайны степени m, которые имеют отличные от нуля значения на (m+1) участках сплайна, называемых интервалом носителем B-сплайна. Значения B-сплайна на границах участков называются его узлами.

В общем случае каждый участок базисного B-сплайна может включать d интервалов Δt=t_i – t_i+1 (d=1, 2, 3, …). При этом величина интервала носителя B-сплайна равна d∙Δt·(m+1), а значение сплайна степени m в момент времени t_i определяется по формуле , (2), где l – число узлов базисного B-сплайна на половине интервала носителя, в которых его значения отличны от нуля; – величина B-сплайна степени m в k-м узле его интервала носителя.

Величина l зависит от степени сплайна и для сплайнов нечётной степени m определяется выражением .

Для определения полиномиальных сплайнов (1) применяются разные способы вычисления коэффициентов C_i [1]. Рассмотрим случай, когда коэффициент C_i, соответствующий i-му В-сплайну , определяется в виде взвешенной суммы произведений отсчётов исходного сигнала, которые попадают в интервал носитель i-го B-сплайна, на соответствующие значения B-сплайна (далее индекс m в обозначении B-сплайна опущен) . (3).

Таким образом, с учётом формулы (2) . (4)

Можно показать, что если коэффициенты базисных B-сплайнов полиномиального сплайна произвольной нечётной степени m, определяются по формуле (3), то его амплитудный спектр в отличие от амплитудного спектра исходного сигнала, имеет затухающий характер

. (5)

Следовательно, такая сплайн-аппроксимация может рассматриваться как низкочастотная фильтрация сигнала.

В формуле (4) для определения i-го значения сплайна используется (k·d+n) значений исходного сигнала, где

Введём индекс j=k∙d+n. С учётом величины l и диапазона изменения индексов k и n диапазон изменения индекса j составит

Таким образом, принимая и учитывая свойство симметрии B-сплайнов, формулу (4) можно представить в виде .

Применяя перестановку операций суммирования, получим , (6)

причём, , если .

Последнее условие следует из того, что значения любого базисного B-сплайна за пределами его интервала носителя тождественно равны нулю.

Недостатком формулы (6) является то, что для определения значения сплайна S^m(t_i), соответствующего моменту времени t_i, требуется знать значение исходного сигнала в момент t_i, а также L значений исходного сигнала до момента t_i и L значений после момента t_i. Формула для определения значения сплайна S^m(t_i) по значениям исходного сигнала и сплайна, предшествующим моменту t_i, может быть получена путём следующих преобразований.

Обозначим в формуле (6) , при , , при .

где , если .

Получим . (7).

Из формулы (7) определим . (8)

. (9)

Умножив обе части выражения (8) на v_{L ,} а обе части выражения (9) на v_L-1, найдём их разность

.

Обозначив , получим

. (10).

В последнем выражении отсутствует слагаемое, использующее значение исходного сигнала в момент t_i+(L-1). Применяя подобные преобразования и вводя коэффициенты , r = 1, 2, …, L-1, где r – номер цикла, можно за (L-1) циклов получить выражение, в котором будут отсутствовать значения исходного сигнала после момента t_i.

Таким образом на (L-1) цикле формула (10) преобразуется к виду

.

Из последнего выражения находим

. (11)

С учётом значений коэффициентов выражение (11) принимает следующий вид

, (12), где - коэффициент в степени j.

Можно показать, что коэффициент имеет малую величину, которая уменьшается с увеличением степени m сплайна.

Введя в выражении (12) обозначения

, для; , для. и пренебрегая последним слагаемым, получим .

Очевидно, что коэффициенты a_j, b_j полностью определяют значение сплайна S^m(t_i) в момент времени t_i по значениям сигнала и сплайна, предшествующим моменту t_i. При этом обеспечивается фильтрация сигнала с амплитудно-частотной характеристикой (5).

Литература

1. 
   Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. – М.: Наука, 1980. – 352 с.
   


ALGORITHM OF DIGITAL SIGNALS'FILTRATION ON THE BASIS OF SPLINE APPVOXIMATION METHODS

Prosochkin A.

The problem of the digital filtration of a signal given in the form of discrete readings against uniform time scale by restoring signal with the help of polynomial splines to a m arbitrary degree is being discussed in this paper. Polynomial splines’ representation in the form of weighed basic B-splines’ sum along with the proper coefficients used for their determination Separate B-splines’ coefficient is calculated as weighed sum of signal reading products into appropriate B-spline values. The formula showing that method of basic B-spline coefficients makes its amplitude spectrum damping, is presented in this paper. Consequently, such spline approximation may be considered like low-frequency signal filtration.

The analysis of the presentation from of polynomial splines shows that in order to determine a spline value corresponding a time moment one should know the value of the initial signal at a moment plus a particular number of the initial signal data values before or after moment. Above mentioned splines’ presentation may be used only for the secondary data file processing but not for the real-time signal processing when the result of one processing cycle precedes the following one.

The transformations that help to proceed to new polynomial splines’ representation form are also presented in this paper. Algorithm of weight coefficients’ calculations has been worked out. Formula that determines spline-value at a time – moment within some values of initial signal and spline values preceding to a – moment has been got. That formula realizes recursive algorithm of low-frequency signal filtration and it meets well-known frequency-amplitude characteristic.



ВЛИЯНИЕ СПЕКТРАЛЬНОГО ПРОСАЧИВАНИЯ НА ПОВЕДЕНИЕ АВТОКОРРЕЛЯЦИОННОЙ ФУНКЦИИ УСЕЧЕННОГО гармонического сигнала

Ханян Г.С.

Центральный институт авиационного моторостроения им. П.И. Баранова

111116,  Москва, ул. Авиамоторная, 2. [email protected]

Корреляционная функция используется во многих областях прикладной математики и физики: в теории вероятностей и математической статистике – для оценки характеристик случайных процессов, в экономике – для прогнозирования промышленного развития, в аэрогидродинамике – для определения интегрального массштаба турбулентности и т.д. Применяется корреляционная функция и в алгоритмах интерполяции данных (пикселов) в двумерной обработке изображений, в PIV-методе оптической визуализации потоков и т.д.

В практике цифровой обработки сигналов дискретную автокорреляционную функцию вычисляют с помощью обратного дискретного преобразования Фурье (1), от квадрата модуля дискретной спектральной функции (2), N-точечной цифровой реализации (3)

аналогового сигнала s(t), наблюдаемого в окне времени длительностью T. При этом дискретизация проводится с частотой F=1/ t=N /T, начиная с некоторого момента времени T₀.

Применение формулы Винера-Хинчина (1) для вычисления корреляционной функции таит в себе две проблемы – методологическую и метрологическую. Первая из этих проблем обусловлена периодичностью цифровой реализации s_n, возникающей при обращении ДПФ (2). Это приводит к разрыву запаздывающей последовательности s_n+k в точке n=N–k, из-за чего прямое вычисление отсчетов C_k через отсчеты цифровой реализации s_n осуществляется по двучленной формуле

, (4), в которой на главную часть, описываемую дискретной функцией R_k, накладывается зеркальная часть R_N–_k, формируемая с правого конца цифровой реализации (3). Этот факт хорошо известен в  литературе (см., например [1] ) как наложение во временной области. Проблема состоит в продолжающейся дискуссии о том, какая из корреляционных функций является «истинной» – круговая C_k или линейная R_k. Более того, открытым остается и вопрос о том, на что делить сумму при определении R_k – на число слагаемых N–k, как в (4), или же на число отсчетов N, как при определении C_k. С вычислительной же точки зрения проблема решена: линейная корреляционная функция выражается через круговую, вычисляемую с применением алгоритма БПФ от дописанной N нулями исходной цифровой реализации сигнала.

Вторая, метрологическая проблема относится к точности представления автокорреляционной функции эталонного сигнала – гармонического колебания, происходящего с амплитудой a₀, частотой f₀ и начальной фазой ₀. Цифровая реализация такого сигнала имеет вид

, (5), где m₀ есть адрес максимального пика в спектре амплитуд сигнала A_m=2|S_m|, m=1,2,…,m₀,…N/2–1; ₀ – дробная добавка к адресу пика; ₀ – начальная фаза сигнала в системе отсчета времени, связанной с окном наблюдения. Параметры эти определяются следующим образом: ; . (6)

В работе [2] показано, что именно параметр ₀, изменяющийся в диапазоне –1/2<₀1/2, является основным фактором, искажающим как спектр амплитуд (занижение максимального пика, сопровождаемое просачиванием боковых составляющих), так и спектр фаз (смещение начальной фазы). Естественно поэтому предположить, что параметр этот может влиять и на структуру корреляционной функции. Целью настоящей работы является исследование характера этого влияния и степени вносимого при этом искажения корреляционной функции сигнала (5).

Для решения поставленной задачи получим аналитические выражения для функций C_k и R_k.

Начнем с линейной корреляционной функции. Подстановка (5) во вторую из формул (4) приводит к следующей зависимости дискретной функции R_k от дискретного запаздывания _k=kt:

. (7).

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

Заменив в (7) k на N–k, приходим к выражению для R_N–k в третьей формуле (4). Подставив обе вычисленные линейные функции в первую из формул (4), получаем точное аналитическое выражение для круговой автокорреляционной функции усеченного гармонического сигнала:

. (8)

В приложениях принято иметь дело с корреляционной функцией, центрированной по сравнению со средним значением и нормированной по отношению к дисперсии цифровой реализации сигнала. В центрированной и нормированной круговой корреляционной функции (9), фигурирует квадрат среднего значения цифровой реализации, равного как раз нулевому отсчету S₀ спектральной функции (2). Вычисление S₀ и подстановка его в первую из формул (9) приводит ее к

, (10), где постоянные (не зависящие от переменной времени _k) константы ₀ и ₀ суть . (11).

Проанализируем полученный результат. Очевидно, что известная в литературе «классическая» формула для автокорреляционной функции гармонического сигнала (12), вытекает из (10) лишь при нулевом значении параметра ₀ – т.е. при отсутствии спектрального просачивания. При ₀0 корреляционная функция искажена, причем характер этого искажения определяется вторым слагаемым в числителе (10), имеющим вид произведения функции от _k на множитель sin ₀, изменяющийся от 0 до 1 при изменении |₀| от 0 до 1/2. Второе слагаемое в знаменателе равно значению этого произведения при _k=0 и есть константа. Константой является и третье, квадратичное по ₀ слагаемое, общее для числителя и знаменателя и представляющее собой квадрат среднего значения S₀, деленный на половину квадрата амплитуды сигнала a₀.

Оценим порядок величин ₀ и ₀. Рассмотрим гармонический сигнал, частота f₀ которого находится в окрестности одной из трех характерных точек частотной оси – в начале спектра ( f₀0), в середине ( f₀F/4) и вблизи частоты Найквиста ( f₀F/2). Пусть f – расстояние частоты от указанной характерной точки, малое по сравнению с частотой дискретизации F.

В первом случае (низкочастотное колебание) принимаем  f₀ = f и для оценки ₀, ₀ имеем

, (13), где m = [ f  T] – бин, соответствующий размерной частоте f , малый по сравнению с N.

Во втором случае (среднечастотное колебание) частота сигнала f₀ = F/4  f может находиться справа или слева от центральной частоты спектра, что приводит к следующей оценке

. (14)

В третьем случае (высокочастотное колебание) в (11) подставляем f₀ = F/2– f и получаем

. (15)

Учитывая, что m<<N и |cos₀|1, |cos2₀|1, мы можем обобщить все три случая (13)-(15), и по самой грубой оценке утверждать, что для всех значений ₀ и ₀ имеют место ограничения

. (16).

На практике обычно N=2048, и для такой длины цифровой реализации |₀|<1%, (₀)²<1% уже при m16 бин, т.е. при отступлении частоты синусоиды от краев спектра всего на 16/10241,6%. Это позволяет пренебречь всеми постоянными слагаемыми в формуле (10) – третьим в числителе и вторым и третьим в знаменателе. Если же провести сравнение _k /T с ₀, то уже при _k /Tk/N>₀ можно пренебречь и влиянием ₀ на поведение числителя. Для этого, при только что указанных условиях должно быть k>N₀=20481%20, т.е. после 20-го отсчета корреляционной функции вместо точного аналитического выражения (10) можно пользоваться приближенной формулой

, (17), показывающей что при наличии спектрального просачивания огибающая корреляционной функции гармонического сигнала меняется по линейному закону с коэффициентом sin ₀.

Тот факт, что эффект искажения корреляционной функции (12) состоит в ее сужении, не совсем очевиден из (17) – из-за того, что оба осциллирующих с частой f₀ слагаемых отличаются по фазе. Сужение наглядно проявляется при численном моделировании формулы (10). Оно тем сильнее, чем ближе значение |₀| к 1/2 и максимально при наибольшем просачивании амплитуды, т.е. при ₀ = 1/2. При любом из этих двух экстремальных значений параметра ₀ формула (17) принимает простой вид , (18), откуда видно, что корреляционная функция сужается от единицы (при _k=0) до нуля (при _k=T/2).

На рис.1 приведены результаты спектрального и корреляционного анализа гармонического сигнала (5), моделируемого при следующих значениях параметров: a₀=100, m₀=16, ₀=/4, T₀=0, N=512, T =1000 мс. Анализ проводился при ₀=0 (верхнее окно) и ₀=1/2 (нижнее окно). Корреляционная функция (9) вычислялась как с применением ДПФ (1), так и моделированием формулы (10). Видимых различий при этом обнаружено не было, поэтому в окнах представлены результаты применения ДПФ. Видно влияние параметра ₀ как на структуру спектра амплитуд (занижение и уширение пика), так и на структуру корреляционной функции (сужение огибающей).

Продемонстрируем теперь влияние шума на эффект сужения круговой корреляционной функции. На рис.2 приведены результаты спектрального анализа чистого и зашумленного гармонического сигнала (с полуразмахом равномерно распределенной шумовой компоненты, равной амплитуде сигнала a₀). Параметры сигнала те же, что и на рис.1, за исключением того, что спектрограммы и коррелограммы получены для ₀=1/3. Видно, что шум и спектральное просачивание влияют на структуру корреляционной функции независимо друг от друга. Первый из этих факторов делает нулевой и начальные отсчеты корреляционной функции доминирующими над остальными, второй – сужает ее огибающую.

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

Рис.1. Корреляционная функция и спектр амплитуд чистого и зашумленного гармонического сигнала

Вернемся к линейной корреляционной функции. Поскольку центрирование, как было показано для круговой функции, не является значимым искажающим фактором, для простоты анализа рассмотрим только нормированную линейную функцию, полученную из (7): .(19)

Максимальное значение множителя перед ₀ достигается при максимальном значении переменной _k=_N/2=T/2 и равно 2. Поэтому, руководствуясь оценкой (16), можно с установленной выше точностью считать, что линейная корреляционная функция практически не подвержена влиянию эффекта просачивания и для гармонического сигнала имеет вид косинусоиды , (20), о чем свидетельствуют практически неотличимые результаты моделирования формул (19) и (20). Ясно теперь, что при делении суммы в (4) на N вместо N–k функция R_k будет сильно искажена.

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

Литература

1. 
   J.S.Bendat, A.G.Piersol, Engineering applications of correlation and spectral analysis, New York, 1980.
   
2. 
   G.S.Khanyan, “Development of methods for spectral estimation of noisy harmonic signal parameters,” in Proc. 7th Int. Workshop on Spectral Methods and Multirate Signal Processing, Moscow, Sep. 2007,  Tampere International Center for Signal Processing, TICSP series 36, pp. 84-87.
   

