Радиотехника и электроника, 2021, T. 66, № 5, стр. 443-456

ВВЕДЕНИЕ

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

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

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

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

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

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

Известные по литературе функции, удовлетворяющие этим требованиям имеют вид начального отрезка ряда Фурье с небольшим количеством членов, что обеспечивает простоту вычислений. Малое количество членов ряда позволяет, тем не менее, получить высокие характеристики, поэтому такие функции наиболее востребованы в приложениях. Возможные значения асимптотической (с ростом частоты) скорости спада лепестков для таких функций ограничены набором чисел 6(2k + 1), k = 0, 1, …. Таковыми функциями являются функции Хемминга [1] (одна функция), Блекмана–Хэрриса [1] (одна функция), Наттелла [2] (три функции). В работах Албрехта [3] и независимо Зайцева [4, 5] разработаны методы синтеза континуального семейства оптимальных по минимаксному критерию функций указанного вида с различным уровнем подавления боковых лепестков и табулированы их параметры.

Наиболее широкое семейство оптимальных функций, удовлетворяющих сформулированным выше требованиям, получено в работах авторов [6–10], где описан метод синтеза и приведены таблицы параметров большого количества оптимальных весовых функций, представимых в виде линейной комбинации малого числа членов, каждый из которых является степенью одной и той же косинусоидальной функции. Разработанный метод позволяет синтезировать оптимальную функцию по следующему набору параметров: число членов ряда, ширина главного лепестка спектра и асимптотическая скорость спада его лепестков V ≥ 6 дБ/окт. Заметим, что в литературе, кроме указанных работ, авторам не известно описания функций со скоростью спада, некратной 6 дБ/окт.

Простейший способ вычисления имеют, очевидно, функции в виде полинома, однако в литературе подобные функции с перечисленными свойствами не описаны. Получено лишь незначительное количество функций [11, 12] с узким диапазоном параметров, не перекрывающим указанный выше динамический диапазон. Кроме того, эти функции не являются оптимальными по описанному выше критерию.

В настоящей работе предложено семейство полиномиальных весовых функций, имеющих минимальный уровень максимального бокового лепестка и произвольную асимптотическую скорость спада V ≥ 6 дБ/окт. Функции имеют вид линейной комбинации степеней одного и того же полинома второго порядка, характеризуются малой вычислительной сложностью и конструируются путем численного решения минимаксной задачи оптимизации. Приведены параметры репрезентативного множества синтезированных функций.

1. ПРЕДСТАВЛЕНИЕ ОКНА И БАЗИСНЫЕ ФУНКЦИИ

Весовая функция w(t), используемая при гармоническом анализе непрерывного сигнала на интервале наблюдения –T/2≤ t ≤ T/2 длительностью T, является финитной четной положительной непрерывной функцией времени (кроме, возможно, краев интервала): w(t) = 0 при |t| > T/2.

Спектр W_c(f) весовой функции w(t) определяется интегральным преобразованием Фурье:

При цифровой обработке сигналов окно представлено N временными отсчетами, взятыми с шагом δ на интервале наблюдения сигнала T = Nδ. Расположим начало координат в центре между крайними отсчетами. Тогда моменты выборок равны t_q = – 0.5 (T – δ) + qδ, q = 0, 1, …, N – 1. Дискретную весовую функцию w(t_q) будем рассматривать как отсчеты непрерывной функции w(t) в точках t_q. Спектр дискретного окна определяется формулой дискретного по времени преобразования Фурье [13]

Спектр (2) периодичен по частоте f с периодом 1/δ, и его достаточно анализировать на периоде f ∈[–0.5/δ, 0.5/δ]. При выборе достаточно малого шага δ спектры (1) и (2) различаются незначительно [13].

Ниже решается задача синтеза весовых функций и приводятся таблицы их параметров для случая цифровой обработки. Один из этих параметров, а именно, асимптотическая (при f → ∞) скорость V спада боковых лепестков спектра, относится, строго говоря, к соответствующей непрерывной функции w(t), определяется как величина снижения огибающей модуля спектра W_c ( f ) при увеличении частоты вдвое и измеряется в децибелах на октаву. Тем не менее, этот параметр оказывается полезен для описания и дискретной весовой функции, когда число выборок N велико, а частота не приближается к половине периода 1/2δ. Для краткости слово “асимптотическая” в дальнейшем опускается.

Рассматриваемые в данной работе весовые функции имеют вид линейной комбинации небольшого числа m + 1 базисных функций ψ(p_k, t):

где b₀, b₁, …, b_m – коэффициенты, а p_k – параметр, определяющий форму k-той базисной функции.

Базисные функции, предлагаемые в настоящей работе, имеют следующий вид:

где p ≥ 0 – неотрицательное действительное число.

Рассмотрим свойства спектра (1) функции (4). Так как функция (4) действительна и четна, то экспоненту в (1) можно заменить косинусом:

Спектр (5) симметричен по частоте относительно начала координат, поэтому ниже без ограничения общности полагаем, что частота положительна: f ≥ 0 .

Интеграл (5) согласно [14, п. 3.387.2], при p > –1 равен

где Г(∙) – гамма-функция, J(∙) – функция Бесселя.

При больших значениях аргумента функция Бесселя аппроксимируется формулой [14, п. 8.451.1]

Поэтому при f → ∞ спектр (6) базисной функции имеет асимптотическое выражение

Отсюда следует, что при f → ∞ спектр (6) имеет асимптотическую оценку O(f^–(p+1)). Следовательно, огибающая боковых лепестков спектра базисной функции (4) спадает со скоростью V = 6(p + + 1) дБ/окт. Выбирая степень p в соответствии с соотношением p = V/6 – 1 можно получить произвольно заданное значение скорости спада V ≥ 6 дБ/окт. Насколько известно авторам, весовые функции со скоростью спада лепестков, некратной 6 дБ/окт, в литературе не описаны.

Нули спектра (6) совпадают с нулями функции Бесселя J_{p + 0.5} (πfT). При нумерации нулей от 1 асимптотическое положение нуля f_zn с большим номером n следует из (7):

Весовая функция (3) содержит (m + 1) базисных функций вида (4). Для получения низких боковых лепестков в ее спектре лепестки одной базисной функции должны компенсироваться лепестками другой. Для хорошей компенсации целесообразно сделать так, чтобы нули спектров последовательных базисных функций совпадали, а их боковые лепестки были противофазны. Из (8) следует, что асимптотически это достигается, если разности степеней p₀, p₁, …, p_m базисных функций, используемых в (3), равны целым четным числам. Это условие удовлетворяется при выборе степеней в соответствии с соотношением p_k = μ + 2k, k = 0, 1, …, m, где μ – неотрицательное действительное число. При этом весовая функция (3) приобретает вид

Величину m будем называть порядком весовой функции. Будем также считать, что общий множитель в этой формуле вынесен за знак суммы, т.е. b₀ ≠ 0. Тогда скорость спада боковых лепестков спектра функции (9) определяется параметром μ и составляет V = 6(μ + 1) дБ/окт.

Степень μ множителя перед суммой в (9) определяет порядок нуля функции на краях ее области определения и скорость убывания лепестков спектра V = 6(μ + 1) дБ/окт, а коэффициенты суммы подбираются ниже так, чтобы уровень максимального лепестка спектра окна был наименьшим.

Сумму в (9) можно рассматривать как полином с четными степенями относительно переменной x = 1 – (2t/T)². Вычисление функции вида (9) требует малого объема вычислений, поскольку число членов в сумме (9) мало, выражение в квадратных скобках одинаково для всех членов, а сумма может быть вычислена по схеме Горнера.

Для целых µ выражение (9) можно записать в традиционной форме в виде полинома относительно переменной (2t/T)². Однако приведенная форма более наглядна и имеет ряд описанных ниже преимуществ.

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

2. МИНИМАКСНАЯ ЗАДАЧА

Рассмотрим метод решения задачи синтеза весовых функций для дискретного случая. Для решения необходимо надлежащим образом выбрать неизвестные коэффициенты функции (9). Объединим эти коэффициенты в вектор b = (b₀, b₁, …, b_m) и включим его в явном виде в список аргументов весовой функции (9) и ее дискретного спектра (2): w(b, t) = w(t) и W(b, t) = W(t).

Подставив выражение для весовой функции (9) в формулу (2) для спектра, получаем

где Ψ(μ + 2k, f ) – дискретный спектр базисной функции ψ(p_k = μ + 2k, t).

Пусть W_max > 0 – максимальное значение модуля спектра (10). Определим функцию качества φ(b) как максимальное значение нормированного модуля спектра |W (b, f)/W_max| вне заданной области частот [–β/T, β/T], где β/T в настоящей работе определяет максимально допустимую полуширину главного лепестка по уровню максимального бокового лепестка:

Величина 1/T называется бином [1], при этом полоса β/T составляет β бин.

Построение оптимальной весовой функции заключается в отыскании вектора коэффициентов b, доставляющего минимум ρ целевой функции φ(b):

где правая часть означает аргумент b функции φ(b), при котором достигается минимум.

В работах авторов [6–10] разработан метод решения подобной задачи для случая тригонометрических базисных функций. Этот метод пригоден фактически без изменений и для рассматриваемых в данной работе базисных функций (4). В связи с этим ограничимся здесь только кратким перечислением основных этапов метода. Подробное описание может быть найдено в [8, 10].

1. Преобразование задачи к линейному виду.

Рассмотрим свойства функции вида (9). Во-первых, максимум дискретного спектра (2) достигается при f = 0. Действительно, в выражении (2) для спектра суммируются неотрицательные числа w(b,t_p) с весом cos(2πt_pf). Очевидно, максимум достигается, если все веса равны единице, что реализуется только при f = 0.

Во-вторых, умножение всех коэффициентов b₀, b₁, …, b_m на фиксированный множитель не изменяет функцию качества (11). Для обеспечения единственности решения зафиксируем максимальное значение спектра W_max = W(b, 0) =1.

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

2. Анализ функции качества.

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

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

3. Выбор алгоритма оптимизации.

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

4. Верификация полученного решения.

Так как используемые алгоритмы поиска оптимума являются численными, то необходимо проверить, насколько полученный результат близок к точному решению. Авторами разработан метод такой проверки [8, 10] на основе использования нижней границы минимума, полученной в [15]. Метод дает возможность получить оценку близости полученного решения к оптимальному.

3. КАТАЛОГ СИНТЕЗИРОВАННЫХ ОПТИМАЛЬНЫХ ВЕСОВЫХ ФУНКЦИЙ

Приведем параметры нескольких подклассов весовых функций с различными скоростями спада лепестков, синтезированных по описанной методике. Эти параметры (табл. 1–5) являются наиболее полезными результатами работы, так как для использования весовой функции необходимы только ее коэффициенты. Для сокращения объема таблиц приняты во внимание следующие обстоятельства.

Во-первых, характеристики оптимальных весовых функций слабо зависят от числа выборок N. Для различных N разница в уровне максимального лепестка составляет доли, максимум единицы, децибел и тем меньше, чем больше сравниваемые значения N. Для N = 1024 разница с асимптотическим случаем N → ∞, как правило, менее 0.1 дБ. Поэтому в таблицах приводятся характеристики оптимальных функций только для одного значения N = 1024, часто используемого в приложениях. При таком значении N, как коэффициенты, так и характеристики весовых функций, практически совпадают с величинами для непрерывных функций.

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

Уменьшение уровня максимального бокового лепестка при синтезе достигается за счет увеличения β и порядка функции. Расчеты показывают, что при постоянном β увеличение порядка функции сверх минимально необходимого для обеспечения отсутствия нулей функции на интервале [0, β] приводит к довольно слабому улучшению характеристик функции. Поэтому в табл. 1–5 содержатся параметры функций только для указанного минимального порядка.

В таблицах даны следующие общепринятые параметры весовых функций [1]: уровень P максимального бокового лепестка относительно главного лепестка; потери Π_w в отношении сигнал/шум, вызванные использованием весовой функции; потери Π_s за счет рельефа спектра (scalloping loss); эквивалентная шумовая полоса L_n, ширины главного лепестка L₃ и L₆ по уровням минус 3 и минус 6 дБ, соответственно. Так как ширина лепестков спектра обратно пропорциональна T, то все полосы частот и параметр β, как это принято в литературе [1], выражены в относительных единицах (в бинах), т.е. умножены на T.

В табл. 1–5 значения уровня максимального лепестка округлены с точностью до 0.1 дБ, а разрядность коэффициентов выбрана такой, что получающиеся значения P отличаются от оптимальных не более, чем на 0.05 дБ, что гарантировано алгоритмом верификации (этап 4 предыдущего раздела).

Наибольший интерес представляют функции с целым значением параметра μ, так как их вычисление требует наименьшего количества операций. В табл. 1–5 даны параметры синтезированных функций со скоростями спада лепестков 6, 12, 18, 24 и 30 дБ/окт., соответственно. Шаг по параметру β для всех таблиц составляет 0.25 бин. Функции с максимальным лепестком более минус 40 дБ не представляют значительного интереса и поэтому для краткости не приведены.

Функции со скоростью спада лепестков, некратной 6 дБ/окт, не имеют аналогов в литературе, кроме построенных авторами косинус-полиномиальных функций [9, 10]. Изложенный выше метод позволяет синтезировать полиномиальные весовые функции и с этими скоростями спада лепестков, однако в данной работе такие функции не приводятся поскольку при сопоставимых характеристиках они требуют большего объема вычислений, чем функции со скоростями кратными 6 дБ/окт.

4. ХАРАКТЕРИСТИКИ ОПТИМАЛЬНЫХ ВЕСОВЫХ ФУНКЦИЙ

Рассмотрим характеристики построенных полиномиальных функций. Рисунок 1 иллюстрирует зависимость уровня максимального бокового лепестка P от параметра синтеза β при скоростях спада боковых лепестков V = 6, 18 и 30 дБ/окт. Эти значения выбраны таким образом, чтобы не произошло слияния кривых на рисунке.

Рис. 1.

Зависимости уровня P максимального бокового лепестка от параметра β для оптимальных функций со скоростями спада лепестков 6 (1), 18 (2) и 30 (3) дБ/окт.

Из рис. 1 видно, что увеличение скорости спада лепестков при постоянном β приводит к некоторому возрастанию уровня максимального бокового лепестка. Для фиксированной скорости спада каждая кривая хорошо аппроксимируется прямой линией, небольшие отклонения от которой имеют случайный вид. Для грубой оценки значения параметра синтеза β, требуемого для получения желаемого уровня P максимального бокового лепестка можно использовать линейное соотношение P ≈ (–27.7β + 19.8) дБ. При увеличении скорости спада лепестков увеличивается крутизна аппроксимирующих прямых с ‒ 27.3 дБ/бин для V = = 6 дБ/окт до ‒ 29.7 дБ/бин для V = 30 дБ/окт. Это приводит к сближению кривых между собой при увеличении β.

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

Рисунок 2 иллюстрирует зависимости шумовой полосы L_n и ширины главного лепестка L₆ по уровню минус 6 дБ от ширины главного лепестка L₃ по уровню половинной мощности. На рисунке приведены указанные параметры более, чем 200 оптимальных функций. Каждая функция отображается одной точкой для каждой из зависимостей.

Рис. 2.

Зависимости шумовой полосы L_n (1) и полосы L₆ (2) от полосы L₃.

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

С использованием этих уравнений можно довольно точно оценить полосы L_n и L₆, по значению полосы L₃. Кроме того, полоса L_n определяет потери в отношении сигнал/шум, обусловленные введением весовой функции: Π_w = 10 lg L_n. Таким образом, из четырех перечисленных параметров три являются зависимыми.

Малая величина отклонения реальных полос оптимальных весовых функций от прямых линий позволяет использовать эти отклонения в качестве критериев качества весовой функции.

Проведем сравнение весовых функций, полученных в данной работе, с известными. Так как параметры весовых функций зависят от скорости спада лепестков, то это сравнение следует проводить для одной и той же скорости спада лепестков. На рис. 3 такое сравнение дано для скорости спада 6 дБ/окт, где показана взаимозависимость двух основных параметров: ширины главного лепестка по уровню половинной мощности и уровня максимального бокового, для нескольких классов классических весовых функций: Кайзера–Бесселя, Гаусса, Блекмана–Хэрриса [1] наряду с новыми функциями, построенными в данной работе.

Рис. 3.

Зависимости ширины L₃ главного лепестка по половинной мощности от уровня P максимального бокового лепестка для функций Гаусса (1), Кайзера–Бесселя (2), Блэкмана–Харриса (3, отдельные кружки ), Дольфа–Чебышева (4) и функций табл. 2 (5).

Также представлены функции Дольфа–Чебышева [1], которые имеют нулевую скоростью спада лепестков, обеспечивают минимально возможную ширину главного лепестка при заданном уровне максимального бокового, поэтому соответствующая им линия является границей достижимых параметров. По объему вычислений функции Дольфа–Чебышева не могут конкурировать с функциями малого порядка, так как имеют сложный алгоритм вычисления и максимальный порядок (примерно N/2), что приводит к несоизмеримо большему объему вычислительных операций по сравнению с функциями данной работы, порядок которых не превосходит 12.

Чем ниже расположена кривая на графике рис. 3, тем выше качество соответствующего класса весовых функций с точки зрения решаемой задачи, т.е. наглядно видно преимущество построенного класса функций над известными. Для других скоростей спада лепестков также наблюдается аналогичное преимущество. При фиксированном уровне бокового лепестка P < –40 дБ, ширина главного лепестка по половинной мощности L₃ новых функций, несмотря на малый порядок, всего на 2…4% превышает граничное значение.

Кривые для полиномиальных функций на рис. 1 и 3 имеют небольшие выбросы вверх. Соответствующие функции не приведены в таблицах и имеют максимально возможный параметр β для данного порядка функции, при котором отсутствуют нули спектра в области [0, β]. Эти функции имеют также минимально возможный для данного порядка уровень максимального лепестка. Выбросы указывают на небольшое ухудшение качества соответствующих оптимальных функций.

Сравнение описываемых в данной работе функций с косинус-полиномиальными функциями, синтезированными авторами в работах [4‒10], показывает, что их характеристики весьма близки при одинаковых параметрах β и μ. Поэтому на рис. 3 кривая для косинус-полиномиальных функций не приведена во избежание наложения на кривую для полиномиальных функций.

Цена за упрощение вычислений для полиномиальных функций состоит в увеличении порядка по сравнению с косинус-полиномиальными функциями с одинаковыми характеристиками. На рис. 4 показана зависимость порядка оптимальной функции от значения максимального бокового лепестка для полиномиальных (1) и косинус-полиномиальных (2) функций при скорости спада лепестков 6 дБ/окт. При значительных боковых лепестках порядки отличаются незначительно, однако при малых уровнях порядок полиномиальных функций примерно в 1.5 раза больше. Такая ситуация имеет место и для других скоростей спада лепестков.

Рис. 4.

Зависимости порядка оптимальной весовой функции от уровня максимального лепестка для полиномиальных (1) и косинус-полиномиальных (2) функций.

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

Сравнение по объему вычислений с известными полиномиальными функциями [11, 12] показывает преимущество построенных в данной работе функций, так как при сопоставимых максимальных уровнях боковых лепестков известные функции имеют более высокий порядок.