CC BY
90
Электронное научно-техническое издание
НАУКА и ОБРАЗОВАНИЕ
Эя №ФС 77 - 30569. Государственная регистрация №0421100025.155Н 1994-0408_
Теоретические основы спектрального анализа в базисе Хартли 77-30569/230816
# 10, октябрь 2011 Сюзев В. В.
УДК 519216.1/2
МГТУ им. Н.Э. Баумана [email protected]
Условные обозначения:
БП - Быстрых преобразований.
БПФ - Быстрых преобразований Фурье
ДЭФ - Дискретные комплексно-экспоненциальных функций;
НИР - Научно-исследовательская работа.
ФЦП - Федеральная целевая программа
Широкая сфера применения комплексно-экспоненциальных базисных функций в различных областях науки и техники общеизвестна [1]. Тем не менее базис этих функций обладает существенным недостатком, связанным с его комплексным характером. Хорошей альтернативой ему служит базис Хартли, образованный из тех же тригонометрических функций, что и комплексно-экспоненциальный базис, но принимающий только вещественные значения [2]. Основам спектрального анализа в базисе Хартли и посвящена данная статья, в которой наряду с известными сведениями приводятся и новые оригинальные результаты теории функций и преобразований Хартли.
1. Функции и преобразования Хартли
В 1942 г. Р.Хартли ^.^г^у) предложил для спектрального анализа систему базисных функций {cas(kz)}, каждая из которых являлась суммой соответствующих тригонометрических функций [2]:
са8(к2) =со$(к2)+$1п(кг).
Эта система представляла собой по сути формальную базисную систему, образованную из четных 2со&(кг) и нечетных 2sin(kz) простых тригонометрических функций, что позволило
в данном случае из ненормированных функций cos(kz) и sin(kz) получить систему нормированных функций (са$^)}.
Базис Хартли является системой вещественных периодических функций, определенных на любом интервале длительности 2п с периодом, так же равным 2п. Он является ортонормированным на том же интервале, т.к. для любых его функций
— í cas(kz)cas(mz)d 2n ¿ I 0, k Ф m,
[1, k = m,
Z = ■
и может быть использован для разложения в ряд Фурье-Хартли математических функций, определенных на интервале, кратном 2п, и удовлетворяющих на нем условию интегрируемости с квадратом.
Для представления непрерывных сигналов с конечным временным интервалом определения длительностью Т базис Хартли должен быть преобразован к виду
f 2п Л f2п1 • f 2п1 c a s I — kt I = cos I — kt 1 + sin I — kt I
l T J l T J 1 T J
или
c a s (2n ft) = cos (2n ft) + sin (2n ft),
где f=k/T является частотой. Пара непрерывных преобразований Хартли записывается следующим образом:
X (k ) = T í x(t )ca s I fkt "jdt,
x(t) = ¿ X(k)cas I ^kt 1,
k=0 l T J
(1)
а соответствующее им равенство Парсеваля имеет следующий вид:
1 . ш
-1 х2 а ^х).
Т Т к=0
Ряд Хартли (1) для сигналов, принадлежащих пространству Ь2р, обладает
среднеквадратической сходимостью к этим сигналам. Спектр Х(к) Хартли вещественных сигналов так же является вещественным.
Все родственные системы: тригонометрическая, комплексно-экспоненциальная и Хартли содержат одни и те же простые функции. Поэтому между спектрами сигнала в этих базисах существует определенная связь. Действительно имеем
XХ(0)=XЧ(0),
Хх (к) = [Х¥ (к) + Хн (к)]/ 2, к Ф 0, (2)
Хх (к ) = Хэ (к )] -1 т [ Хэ (к )].
Здесь в виде Хх(к) обозначен спектр Хартли, в виде Хэ(к) - спектр в комплексно-экспоненциальном базисе, а для тригонометрического спектра используется его обозначение в виде ХЧк) и Хн(к) [3].
Обратная связь спектров базируется на свойстве симметрии тригонометрических функций. Если представить функцию ХХ(к) в виде суммы четной и нечетной компонент Ххч(к) и ХХН(к) соответственно, то
(к) = [ Хх (к)+Хх (- к)] /2 = Т \ ) кг
Ххн (к) = [Хх (к) - Хх (-к)] / 2 = Т \ х() кг)йг.
Тогда будут справедливы следующие соотношения:
Хч (к) = 2(к), Хн (к) = 2Ххн (к),
(3)
Re[ Хэ (к)] = Ххч (к), М Хэ (к)] = - Ххн (к).
Соотношения (2) и (3) взаимосвязи спектров в родственных базисах могут оказаться полезными при спектральном анализе. Например, они позволяют выразить энергетический и фазовый спектры сигнала в комплексно-экспоненциальном базисе через спектр Хартли:
Re2 [Хэ (к)] + 1т2 [Хэ (к)] = Х2н (к) + Х2Ш (к) = = (1 / 4)[Хх (к) + ХХ (-к)]2 + (1 / 4)[ХХ (к) - ХХ (-к)]2 = = [ Х2х(к) + X2 х (-к)]/2.
^ [Хэ (к)] = аг^[-Ххн (к) / Ххч (к)] = = ат<Л%{ [Хх (-к) - Хх (к)] /[Хх (к) + Хх (-к)]}.
Используя понятие частоты, от непрерывных преобразований Хартли нетрудно перейти к интегральным преобразованиям:
ад
х(/) = | х(г)c а s(2п/г
—ад
ад
х(г) =\ х(/)#,
—ад
которые могут быть использованы при представлении непрерывных сигналов, определенных на бесконечных интервалах времени. Равенство Парсеваля в этом случае приобретает также полностью интегральный вид:
J x\t)dt = J X \f )df.
Спектральная плотность X(f) в этих выражениях является вещественной функцией частоты.
Для интегральных преобразований так же справедливы все соотношения взаимосвязи спектров в родственных тригонометрических базисах.
Спектральное представление дискретных финитных сигналов возможно только по
системе дискретных ортонормированных базисных функций Хартли \ca s | 2nki |L
l lN Л
которые могут быть получены путем дискретизации соответствующих непрерывных функций. Процедура дискретизации в этом случае та же, что используется при дискретизации тригонометрической и комплексно-экспоненциальной систем [3], а условие ортонормированности дискретных функций Хартли на интервале [0,N) выглядит так:
1 N_1 Í2п,Л (2п Л Í1, k = m — > c a s I — ki | c a s I — mi | = \ N i~0 ^ N J ^ N J [O, k * m.
Дискретные преобразования Хартли на этом интервале имеют следующий вид:
X(k) = — V x(i)ca s | — ki |,
x(i) = ^ X(k)cas^— ki |,
(4)
где первое выражение является прямым преобразованием Хартли, а второе - обратным. Равенство Парсеваля
N х х 2(/)=х х \к)
М 1=0 к=0
является математическим и физическим подтверждением эквивалентности временной и спектральной областей представления сигналов и полноты дискретной системы Хартли. Условие ортонормированности функций Хартли, пара дискретных преобразований Хартли и равенство Парсеваля на интервале [-N/2, N/2) имеют аналогичный вид и отличаются от их записи на интервале [0, Ы) только пределами суммирования.
Соотношения взаимосвязи спектров в дискретном варианте также справедливы, причем вид их записи полностью совпадает с соответствующими выражениями (2) и (3). Справедлива также запись энергетического и фазового спектров для дискретных комплексно-экспоненциальных функций (ДЭФ) с помощью дискретных спектров Хартли.
Для функций Хартли не выполняется свойство мультипликативности, поэтому для них в прямом виде теоремы спектров мультипликативных базисов [4] не справедливы. Однако, благодаря связи спектров сигналов в базисах ДЭФ и Хартли, формулировка этих теорем, справедливая для комплексно-экспоненциального базиса, может быть записана с использованием спектров Хартли или, как еще говорят, в терминах спектров Хартли. Покажем это на примере двух теорем - теоремы о сдвиге и теоремы о свертке.
В соответствии с теоремой о сдвиге в базисе ДЭФ частотный спектр Уэ(к) сдвинутого сигнала х(1-10) равен модуляции аналогичного спектра Хэ(к) несдвинутого
сигнала х(г) базисной функцией ехр^, т.е.
¥Э (к) = ХЭ (к) ехР к10 ] •
Представим этот спектр в развернутом виде:
2ж 2ж
^э (к) = {М Хэ (к) + 71т[ Хэ (к)]} • (сов(^- Ч ) - 7 яп^ Ч )) =
2 ж 2 ж
= {Яе[Х э (к)] со8(— к10) + 1т[Х э (к)] яп(— Ц,)} +
2п 2п
+7{1т[Хэ (к)] собС-^ Ч) - МХэ (к)] БШ^ Ч )}.
Теперь учтем связь частотного спектра со спектром Хартли (см. (3)). Тогда
Гэ (к) = Re[Y3 (к)] + j Im[Y3 (к)] = [ (к) cos(N ki0) - (к) sm(N kic)] + +j[-Ххч (к) sin(2n kic ) - ХХН (к) COs(N kic )].
Но
2П
7Х (к) = Re[Y, (к)] - 1т[7э (к)] = [ (к) + Хш (к)] • cos^ Ы{)) +
2п
+[ (к) - Хж (к)] sin^^o) =
= Хх (к) cos(2n ко) + Хх (-к) ко)-
N N (5)
Из этого выражения следует, что при сдвиге сигнала во времени спектр несдвинутого
сигнала модулируется косинусной составляющей функции Хартли cos ^j, а его
зеркальная копия - синусной составляющей функции Хартли sin j.
Результирующий спектр Хартли сдвинутого сигнала получатся простым суммированием модулированных составляющих спектра несдвинутого сигнала. В этом состоит суть теоремы о сдвиге в базисе Хартли.
Теперь перейдем к теореме о свертке. В соответствии с ней в базисе ДЭФ спектр Гэ (к) циклической свертки двух сигналов со спектрами Хэ (к) и иэ (к) с точностью до постоянного множителя N равен произведению этих спектров:
Уэ (к) = NХэ (к) иэ (к).
Развернем эти запись, выделив действительную и мнимую части, и учтем отмеченную ранее связь спектров ДЭФ и Хартли. Тогда после несложных преобразований получим:
Уэ (к) = N[((к)Ux4 (к) - Ххн (к)иш (к)) - j(Х^ (к)иш (к) + ХХН (к)UX4 (к))]
Так как
Х^ (к) = ХХ (к) + ХХ (-к) ; Хж (к) = ХХ (к) - ХХ (-к) ;
_их (к) + их (-к) пл_их (к) - их (-к)
иЧ(к) = ХЧ/ у; иш(к) = ■
то спектр ДЭФ свертки может быть преобразован к следующему виду:
Ys (k ) = N[ Xx (k )UX (- k) + Xx (-k )UX (k)] - jN [(Xx (k )Ux (k) - Xx (-k )Ux (-k )]. По этому выражению легко определяется искомый спектр Хартли дискретной свертки:
Yx (k) = N [Xx (k)Ux (k) - Xx (-k)Ux (-k) + Xx (k)Ux (-k) + Xx (-k)Ux (k)].
2 (6)
Зависимость (6) представляет собой математическое описание теоремы о циклической дискретной свертке в терминах спектра Хартли. Как видно из нее, оно значительно сложнее аналогичного описания в частотной области. И только в частном случае, когда один или оба сигнала, используемые в свертке, являются либо четными, либо нечетными, записи теоремы о свертке в базисах ДЭФ и Хартли совпадают. Это связано с тем, что в этом случае спектры Хартли также являются либо четными, либо нечетными функциями. Так, например, для четного сигнала x(i), для которого x(i)=x(N-i), спектр
N-
2 2 2п
XX (k) = — V x(i) cos —ki
и X(-k)=X(k). Поэтому
Yx (k) = N • Xx (k )[Ux (k) - Ux (-k) + Ux (-k) + Ux (k)] = NXx (k )Ux (k).
Однозначная взаимосвязь спектров в базисах ДЭФ и Хартли позволяет все свойства, присущие частотным спектрам, формулировать и для спектров Хартли. Можно сказать, что нет таких задач, для которых справедливо использование комплексных преобразований Фурье и одновременно не может быть применено вещественное преобразование Хартли. Тем не менее вплоть до настоящего времени область практического применения преобразований Хартли остается существенно более узкой по сравнению с областью приложений преобразований Фурье в тригонометрическом и комплексно-экспоненциальном базисах.
2. Быстрые преобразования Хартли на статических интервалах времени Поскольку базисы дискретных функций Хартли и ДЭФ являются родственными базисными системами, то для функций Хартли, несмотря на их немультипликативность, существуют различные типы быстрых преобразований (БП), подобные известным типам БП Фурье (БПФ) [2,5]. Проиллюстрируем эффективность БП Хартли на примере быстрых преобразований по основанию 2. Алгоритмы БП по другим основаниям могут быть получены аналогичным образом. С ними можно ознакомиться также в работе [5].
Для синтеза алгоритмов БП Хартли потребуются свойства функции Хартли суммы двух аргументов. Запишем эти свойства, используя известные тригонометрические тождества:
c a s(a + в) = cos(a) + sin(a) = cos(a) c a s(P) + sin a c a s(-P), (7)
cas(a) при P = 2mn, m = 0,1,2,...,
cas(a + P) = ■
-cas(a) при P = (2m + 1)n, m = 0,1,2,.... (g)
Перепишем теперь прямое дискретное преобразование Хартли (4) без нормирующего множителя 1/N (N = 2n, n = 1, 2, ...)
N-1
х (к) = V x(i)cas(— ki), к = 0,1, ..., N -1
i=0 N (9)
и применим к нему известные способы прореживания сигнала [4].
Прореживание по времени. Выборка сигнала x(i) на первом уровне прореживания
разделяется на две промежуточные выборки x0 (i) = x(2/) и x1 (i1) = x(2 i1 +1). Поэтому
N/2-1 2П N/2-1 2П
X(к) = V xb(i)cas^к2/1) + V x!(¿l)cas[^— к(2/1 +1)]
i1=0 N i1=0 N
или с учетом равенства (7)
N/2-1 2П 2П N/2-1 2П
X(к) = V x0(i)cas(-ЫЛ + cas(— к) V x1(i1)cos(-kil) +
W Р> 0 N/2 1 N р0> 11 N/2 1
2П N/2-1 2П
+cas(-^7к) V x1(i1)sin^^^kil).
N i=0 N /2 (10)
Вычислим с помощью этого уравнения первую часть спектра X(к1), к1 = 0, 1, ..., N/2-1. Она будет равна
N /2-1 2П 2П N/2-1 2П
Х(к ) = V x0(i )cas(-к i ) + cas(— к ) V x (i )cos(-к i ) +
\1) ^ VN/2 ^ N 1 p 11 N/2 1
2n N/2-1 2П +cas(-V x^^nt^^
Но, поскольку
N /2-1
V x0(i)cas(к) = Х(0)(к),
N/2-1 2П
У х(// )са—— к,/,) = Х(1)(к) £0 N/2 117 4 ^
(Х(0)(к1), Х(1)(к1) - спектры Хартли четной и нечетной промежуточной выборки
соответственно), а
N/2-1 N
У Х(//)со8(-^-*А) = [Х1^) + Х(1)(-к1)]/2 = [Х(Г)(кг) + Х(1)(^-к1)]/2, /1 =0 N' 2 2
(11)
N/2-1 N
У Х(//)8т(-^-к1/1) = [Х^)-Х(1)(-к1)]/2 = [Х1^)-Х1^-к1)]/2, /1 =0 N'2 2
то после преобразований получаем
Хк) = Х(0)(к1) + со8(2Лк1)Х(1)(к1) + 81п(2лк1)Х(1)(N - к1).
Для второй половины спектра Х(к1 + N/2), к1 = 0, 1, ..., N/2-1 из (10) имеем: Х(к1 + N / 2) = У х0 (//) с а -к1/1 + 2л//) +
=0 N / 2
2П N/2-1 2П + саБ(— к +п) У х1 (1 )соб(-к/ + 2л1) +
N 1 ' р0 11 N/2 11 1
2П N/2-1 2П
+ саБ[-(— к +л)1 У х (/1Ып(-кI + 2П ).
N 1 ^ N/2 11 1
Учитывая в этом уравнении свойство функций Хартли (8), периодичность тригонометрических функций синуса и косинуса и соотношения (11), после преобразований получим
Х(к1 + N /2) = Х(0) (к1) - сов(2пк1)Х(1) (к1) - яп^к1)Х(1) (N - к1).
Объединяя уравнения для первой и второй части спектра, приходим к общему алгоритму БП Хартли по основанию 2 на первом уровне прореживания по времени:
Х (к1) = Х (0)(к1) + [сов^) Х (1)(к1) + яп(^Л) Х (1)( N/2 - к1)],
Х (к1 + N /2) = Х(0) (к1) - [соб(— к1) Х(1) (к1) + яЦ— к1) Х (1)( N /2 - к1)],
^^ N (12)
к1 = 0,1, ..., N/2-1
В этом алгоритме спектр исходной выборки сигнала выражается через спектры промежуточных выборок.
Поскольку N /2 делится на 2, то полученный алгоритм можно применить и для вычисления промежуточных спектров, введя новый уровень прореживания. На т-м уровне алгоритм БП Хартли будет иметь следующий вид
X (^..^(кт ) = X ^..^,0)(кт ) + [С08(—т+1 кт ) X ^ ^-Д,,-1 Л^ ) +
2п
+ 8Ш(--т+Г кт )X^..Л^1^ - кт )],
2Пк
X (Я1,Я2,...,Ят-1)(кт + 2П-т ) = X (^2,-,^-1,0)(кт ) - [С08(^-^) X (^!.^1.-.^т-1,1)(кт ) +
2пк
+ X (^2,,^,1)(2"-т - кт )],
кт = 0, 1, ..., 2"-т -1; 0,1; а = 1,2,..., т-1,
(13)
и подобен алгоритму Кули-Тьюки с прореживанием по времен [4]. На последнем (п -1)-м уровне прореживания при полном алгоритме БП Хартли 2"-т+1 = 4, кп-1 = 0,1 и тригонометрические множители принимают только значения 1 или 0. Поэтому при т = (п -1) алгоритм (13) упрощается
X(^1.^2. - ^»-2)(0) = X^•^.■■■^л-2,0)(0) + X(Х1,Х2,.А-2,1) (0) ^^(1) = X(^1'^2.---^„-2,0)(1) + X(^1'^2.---^»-2,1)(1) X(^1.^2. - ^»-2)(2) = X 1^1^2.-^»-2,0)(0) - X(Х1'Х2.---Х"-2,1) (0) X(\Д2.--Л-2)(3) = X(^1'^2.---^»-2,0)(1) - X|^1'^2.--Л-2,1)(1)
Спектры x(ЛlД2' Л"-l)(k"-1) промежуточных выборок в этом случае вычисляются с
помощью 2-точечных дискретных преобразований Хартли, которые также имеют простейший вид:
X (^2, .Л-О (0) = х(2"-2 Хп_1 + ••• + 2Х2 + Х1) + х(2"-1 + 2"- 2 Хп_1 + ••• + 2Х2 + \ ),
(14)
X2п-1)(1) = х(2п-2 Хп-1 + • • • + 2+ ) - х(2"-1 + 2"-2 Хп-1 + ••• + 2Х2 + ).
Здесь учтено, что cas(0) = 1, а cas(n) = -1. Эти уравнения задают начальные условия для итерационного процесса вычисления спектра Хартли по БА (13) при m = n -1, n - 2, ..., 1.
Оценим вычислительную сложность алгоритма БП Хартли. Подсчитаем сначала число умножений MБ. На последнем шаге алгоритма умножений нет. На каждом m-м
шаге из всех остальных шагов выполняется N общих умножений, 2m из которых являются тривиальными, т.к. соответствуют нулевым и единичным значениям тригонометрических констант при индексах im и km, равных 0. Поэтому
МБ =£( N - 2m ) = N (n - 2,5) + 2 = N (log2 N - 2,5) + 2.
m=1
Теперь оценим число сложений АБ. Сложения выполняются на всех шагах, причем на каждом m-м шаге из n - 2 первых шагов выполняется по 1,5N - 2m сложений, а на последнем шаге - 2,5N сложений. Общее число сложений будет равно
Аб =£(1,5N - 2m ) + 2,5N = N(1,5n -1) + 2 = N(1,5log2 N -1) + 2.
m=1
Реальное число вычислительных операций можно получить еще меньше, если учесть тривиальные множители, равные 0 и ±1, имеющие место и при ненулевых значениях индексов im и km . Более близкими к реальным являются оценки
МБ = N(n - 3) + 4 = N(log2 N - 3) + 4,
Аб = 1,5N(n -1) + 2 = 1,5N(log2 N -1) + 2.
Следует отметить, что все приведенные оценки МБ и Аб справедливы для
вещественного входного сигнала. В случае комплексного входного сигнала они удваиваются.
Полный алгоритм БП Хартли полезно представлять в виде ориентированного сигнального графа. По сравнению с БПФ сигнальный граф БП Хартли будет содержать дополнительные узлы, зато умножения будут выполняться только на вещественные константы.
Пример 1. Записать полный алгоритм БП Хартли и построить его сигнальный граф для N=8.
Решение. Алгоритм будет иметь два уровня прореживания. По формулам (13), (14) получаем:
- на втором уровне (m = 2; ^ = 0,1; X2 = 0,1; i2, k2 = 0,1) :
^(ь) = {х (0), х (4)}; Х0Д02) = {х (2), х (6)};
х1,0(г2) ={х (1) х (5)}; х1,1 (г2) ={х (3), х (7)};
х (0'0) (0) = х(0) + х(4); х (0'0) (1) = х(0) — х(4);
X(од) (0) = х(2) + х(6); X(од) (1) = х(2) — х(6);
X (1,0)(0) = х(1) + х(5); X (10)(1) = х(1) — х(5);
X (1,1)(0) = х(3) + х(7); X (1,1)(1) = х(3) — х(7);
X (0)(0) = X (0'0)(0) + X (од)(0); X (0)(1) = X (0,0)(1) + X (од)(1);
X(0) (2) = X (0,0) (0) — X(од) (0); X(0) (3) = X(0'0) (1) — X(од) (1);
X(1) (0) = X (1'0) (0) + X(1Д) (0); X(1) (1) = X (1'0) (1) + X(1Д) (1);
X(1) (2) = X (1'0) (0) — X(1Д) (0); X(1) (3) = X (1'0) (1) — X(1Д) (1);
(т = 1; \ = 0,1; к = 0,1,2,3):
на первом уровне
л/2,
X (0) = X(0) (0) + X(1) (0); X (1) = X(0) (1) + [ X(1) (1) + X(1) (3)];
Р2
X (2) = X(0) (2) + X(1) (2); X (3) = X(0) (3) — [ X(1) (1) — X(1) (3)];
X (4) = X(0) (0) — X(1) (0); X (5) = X(0) (1) — [ X(1) (1) + X(1) (3)];
/2"
X (6) = X(0) (2) — X(1) (2); X (7) = X(0) (3) + [ X(1) (1) — X(1) (3)].
Сигнальный граф этого алгоритма приведен на рис. 1. Его структура похожа на структуру графа алгоритма БПФ Кули-Тьюки [4], хотя и содержит два дополнительных узла. На реализацию этого графа нужно затратить 26 сложений и 2 умножения на вещественные константы.
Проверка. Для проверки алгоритма вычислим с его помощью пятый спектральный
коэффициент: X (5) = X (00)(1) + X (01)(1) — — [ X (10)(1) + X (11)(1)] —
72.
Л«>0)-XО^)] = г(0)-х(4) + х(2)-х(6)- XС1)-X(5)]X(3)-X(7)]-
х(1) - х(5)]х(3) - х(7)] = х(0) -42х(1) + X(2) - X (4)+л/2X(5) - х(6). Расчет по прямому алгоритму дает следующий результат:
7
X (5) = У х(/)cas(5ra /4) =
1—0 х(0) + X (1)cas(5л/4) + х(2)cas(5л/2) + х(3)cas(15л/4) +
+х(4)cas(5л) + х(5)cas(25л /4) + х(6)^(30л /4) + X (7)cas(35л /4) = X (0) -Лх (1) + X (2) -
- X (4) + л/2х(5) - X (6)
Результаты совпадают.
Рис. 1. Сигнальный граф полного БП Хартли по основанию 2 с прореживанием по времени для N=8
Прореживание по частоте. Из выборки сигнала на первом уровне прореживания образуются две промежуточные выборки {х^)} и {х(/1 + N/2)}, а затем с помощью уравнения (9) раздельно записываются спектральные составляющие с четными и нечетными номерами. Для k — 2^ получаем
N/2-1 2П N/2-1 2П
X(2k1) = X х(i)cas(N72**) + X х(\ + N/2)cas(+ 2^) =
—0 — 0
N/2-1 2П М /2-1 2П
= Е [х(/) + х(А + N / 2)]с а 8(N/2 Л1г1) = ^ х(0) (/)с а к/
—0 —0
где
х(0) (/) — х(/) + х(/1 + N /2), (16)
а для к — 2к1 +1 с учетом свойств функций Хартли (7) и (8) -
N/2-1 2п 2п X(2к +1) — Е [х(/1) - х(/1 + N /2)]са в(— /) сов(—— А^^!) +
/, — 0
^ 1 N/2
N/2-1 2П 2П
+ Е [х(г1) - х(г1 + N / 2)] сав(- г1) к1г1).
/1—о ^ ^ '2
Используя соотношения (11), последнее уравнение можно преобразовать к виду
N/2-1 2П 2П
X(2к1 +1) — Е [х(г1) - х(/1 + N /2)] сов(— /1)са -- к1/1) +
/1—о ^ ^ '2
N/2-1 2П 2П
+ Е [х(/1)-х(/1 + N/2)]мп(— /1)са8[—— kl(N/2-/1)]
/ —0 ^ ^ ' 2
?
а, изменив с помощью линейной подстановки /1 — N / 2 - /1 порядок суммирования во
второй сумме этого выражения, его можно записать как
N/2-1 2п
Х(2к +1) — Е х(1)(//)са8(-^-к/),
4—0 Я/2 (17)
где
2 п 2п
х(1) (/) — [х(/1) - х(/1 + N /2)] соб(— /1) + [х(N /2 - /1) - х(N - /1)] — /1)
N N (18)
Таким образом, в виде уравнений (15)-(18) получен алгоритм БП Хартли по
основанию 2 на первом уровне прореживания по частоте. Из него следует, что полный
спектр Хартли можно вычислить с помощью N /2 -точечных ДПФ Хартли над суммой и
разностями отсчетов промежуточных выборок, причем разности отсчетов должны быть
предварительно умножены на соответствующие тригонометрические множители.
По способу организации вычислительного процесса этот алгоритм подобен
алгоритму БПФ Кули-Тьюки с прореживанием по частоте [4]. Его можно применить и для
вычисления N /2 -точечных ДПФ Хартли, используя следующий уровень прореживания.
Процедуру прореживания можно продолжать до тех пор, пока число отсчетов в
промежуточных выборках не станет равным 2. В этом случае будет достигнута
77-30569/230816, №10 октябрь 2011 г. Ьпр:/ЛесЬпотац.еёи.ги
14
максимальная глубина прореживания, равная п -1, и получен полный алгоритм БП Хартли с прореживанием по частоте.
На т -м уровне прореживания по частоте БП Хартли по аналогии с БПФ представляется следующим образом:
Х(2т*т + 2т-2 дт_х + ■■■ + 2д2 + дх) = ^ )са.(^т К'т ),
т=о 2
Х(2т£т + 2т-1 + 2т-2 дт-1 +■■■ + 2д2 + д,) = ^ )са8(^т кпгп ), (19)
'т =0 2
кт = 0,1, ..., 2п-т -1; да = 0,1; а = 1,2,..., т -1,
где
X(д1, д2 ,.,дт-1,0) (' ) = х(д1,д2,.,дт-1)(/ ) + ¿(д1 ,д2 ,.,дт-1 ) (' + 2"-т )
\т/ \ т / \ т
( ) 2п
х(«1,«2,...,«т-1,1)(' ) = [х(д1,д2, ,дт-1)(' ) - х ', 2, , т-1 От + 2 )]С08(2П-т+1 'т) +
+[х(д1,д2,..,дт-1)(2п-т - гт) - х (2 'т )]8Ш(2п-т+1) (20)
и описывает итерационный процесс с т = 1, 2,..., п -1 при начальных значениях в виде (16) и (18). На последнем шаге (при т = п -1) ДПФ Хартли в (19) будут 2-точечными и выполняются без умножений. Не потребуются умножения и при вычислении величин (20), поскольку тригонометрические множители в этом случае будут равны 1 либо 0.
Алгоритм БП Хартли с прореживанием по частоте будет содержать такое же число операций, как и алгоритм БП Хартли с прореживанием по времени. Поэтому для него так же справедливы все приведенные ранее оценки. Его можно проиллюстрировать с помощью сигнального графа, который в этом случае будет иметь структуру, близкую к структуре графа БПФ Кули-Тьюки с прореживанием по частоте, хотя и будет содержать дополнительные вычислительные узлы.
Пример 2. Записать алгоритм БП Хартли с прореживанием по частоте и построить его сигнальный граф для N=8.
Решение. Полный алгоритм БП Хартли в этом случае будет иметь два уровня прореживания, т.е т=1,2. В соответствии с общим алгоритмом (19), (20) получаем: - на первом уровне прореживания (т = 1; ' = 0,1, 2, 3; = 0,1):
х(0) (0) = х(0) + х(4); х(0) (1) = х(1) + х(5);
х(0) (2) = х(2) + х(6); х(и) (3) = х(3) + х(7);
х(1) (0) = х(0) - х(4); х(1) (1) = [ х(1) - х(5)] + [ х(3) - х(7)]};
/2
х(1) (2) = х (2) - х(6); х (1)(3) = [ х(1) - х(5)] - [ х(3) - х(7)]};
- на втором уровне прореживания (т = 2; /2, к2 = 0,1; q2 = 0,1):
х(0'0) (0) = х(0) (0) + х(0) (2); х (0'0) (1) = х(0) (1) + х(0) (3);
х (1'0) (0) = х(1) (0) + х(1) (2); х (1'0) (1) = х(1) (1) + х(1) (3);
х (од)(0) = х (0)(0) - х (0)(2); х (од)(1) = х (0)(1) - х (0)(3); х(1Д)(0) = х(1)(0) - х(1)(2); х(1Д)(1) = х(1)(1) - х(1)(3) и результирующий спектр равен:
X (0) = х(0'0) (0) + х(0'0) (1); X (4) = х (0'0) (0) - х (0'0) (1);
X (2) = х(од) (0) + х(од) (1); X (6) = х(од) (0) - х(од) (1);
X (1) = х(1'0)(0) + х (1'0)(1); X (5) = х °'0)(0) - х(1'0)(1);
X (3) = х(1Д) (0) + х(1Д) (1); X (7) = х(и) (0) - х(и) (1).
Рис. 2. Сигнальный граф полного БП Хартли по основанию 2 с прореживанием по частоте для N=8
Сигнальный граф этого алгоритма представлен на рис. 2. Он имеет такие же сложностные характеристики, что и граф рис. 1.
Проверка. Её выполним для того же пятого коэффициента:
X (5) = х (1-0) (0) - х (1-0) (1) = х(1) (0) + х(1) (2) - х(1) (1) - х(1)(3) = х (0) - х (4) + х (2) - х (6) - Т х(1) +
42 42 42 42 42 42 42
х <5) х (3) + ^ х(7) х(з) х (7) х(1) + ^ х(5) = х (0) (1) + х (2) -
- х (4) + л/2х (5) - х (6).
Полученное значение совпадает со значением этого коэффициента, рассчитанным в предыдущем примере.
Используя общую методику синтеза БПФ различного типа и учитывая специфику функций Хартли, аналогичным образом можно записать алгоритмы БП Хартли по другим основаниям, для взаимно-простых множителей и гнездового типа (типа Винограда) [5]. БП Хартли по своим вычислительным характеристикам не уступают БПФ, что только усиливает перспективность их применения для частотного анализа вещественных сигналов.
3. Преобразования Хартли на скользящих интервалах времени
Статические БП Хартли очень близки по описанию к статическим БПФ, что объясняется, как уже отмечалось, родственностью систем функций Хартли и ДЭФ. Это позволяет существенно упростить процедуру синтеза скользящих БП Хартли за счет использования уже разработанных скользящих БПФ. Применим такой подход к синтезу скользящих БП Хартли по основанию 2 с прореживанием по времени.
Сравнение статических БПФ Кули-Тьюки [4] и БП Хартли (13) показывает, что алгоритм Хартли можно получить из алгоритма Кули-Тьюки, если в последнем на каждой итерации одно умножение промежуточного спектра нечетной подвыборки на комплексную константу W-¿и заменить на два умножения одной половины этого спектра
на вещественную константу со8(2лкт / 2" т ), а другой - на вещественную константу 8т(2лкт /2"-т+1). В остальном процедура вычислений остается прежней. Очевидно, что то
же самое можно проделать и применительно к скользящему БПФ Кули-Тьюки с прореживанием по времени [6], получив в результате следующий алгоритм скользящего БП Хартли по основанию 2 с ц уровнями прореживания по времени:
X]m—1\km ) = Xf)(km ) + [^(Пкт / 2^^ ) Xjm¡m-1 к ) +
+ зт(2пкт / 2П—т+1)X(m2m-1(2n—т — кт)],
] 2
^т—1) (кт + 2"—т ) = ^(т) (кт ) — [С0Э(2лкт / 2^'"+1 X',^ (кт ) +
+ зт(2лкт /2П—т+1)X(т>(2п—т — кт)],
3 2
^п—т+1^ у(т) {гуп—т
кт = 0,1, ..., 2п—т —1; т = 1,2,..., д;
XT(h) = х(3 — )со5(2якц1ц /2п—ц).
(21)
При д = п — 1 скользящий алгоритм (21) становится полным.
Для полного алгоритма на последнем (п — 1) -м уровне прореживания тригонометрические константы принимают только значения 1 и 0, поэтому при т = п — 1 алгоритм (21) упрощается
Xf 2)(0) = XJ(n—1)(0) + X——П—2 (0), X<п—2)(1) = хJ(n—1)(1) + X5 п—2П—2(1), XJ(n—2) (2) = X(п—1)(0) — X п ¡1—2(0), X—2)(3)=X(п—1)(1)—X<52(1).
Спектры X3 1)(кп—1) в этом случае вычисляются с помощью 2-точечных дискретных
преобразований Хартли:
X (п—1)(0) = х(3) + х(3 — 2п—1),
X3п—1)(1) = х( 3) — х( 3 — 2п—1).
Эти уравнения определяют начальные данные для итерационного процесса вычисления скользящего спектра Хартли по БА (21) при т = п — 1, п — 2, ..., 1. Его результатом будет
спектр X('p(k), являющийся искомым спектром Xj (к) входного сигнала к з -му моменту
текущему времени.
Проведем оценку сложности скользящего алгоритма (21). На последней итерации умножения не выполняются, а выполняется только 6 сложений, 2 из которых используются при определении начальных данных. На каждой т -й из всех последующих итераций алгоритма выполняются 2 • 2п—т общих умножений и 3 • 2п—т сложений. Однако для двух значений индекса кт, равных 0 и 2п—т—1, тригонометрические функции косинуса
и синуса будут равны 1 или 0 и поэтому в этом случае умножения станут тривиальными, а число слагаемых в уравнениях (21) сократиться до двух. По этой причине число нетривиальных операций на т -й итерации будет равно 2(2п-т - 2) умножений и (3 _ 2п-т - 2) сложений. Отсюда следует и общее число нетривиальных операций в самом алгоритме:
Mc =£ 2(2n-m - 2) = 2(N - 2n) = 2(N - 2 log2 N),
m=1
n-2
Ac (3 • 2n-m - 2) + 6 = 3N - 2(n +1) = 3N - 2(log2 N +1).
m=l (22)
Алгоритм еще можно улучшить за счет дополнительного уменьшения числа умножений, если учесть, что тригонометрические множители при km = 2n-m-2 и
km = 3 • 2n-m-2 по модулю равны между собой и их можно вынести за скобки, обеспечив тем самым выполнение двух умножений вместо четырех. В этом случае на каждой m -й итерации из (n - 2)-х в алгоритме (21) будет выполняться (2n-m+1 - 6) умножений, а их общее число будет равно
Mc = 2[N - (3n - 2)] = 2[N - (3 log2 N - 2)]. (23)
Число дополнительных данных, которые необходимо хранить в памяти для реализации скользящего преобразования Хартли, такое же, как в скользящем алгоритме БПФ Кули-Тьюки по основанию 2 с прореживанием по времени. Пример 3. Записать полный алгоритм скользящего БП Хартли для N=8.
Решение. Алгоритм (21) при n=3 будет иметь всего 2 итерации. Итерация 1 (m=1, k=0, 1, 2, 3):
(0) = *;')(°)+х;«(0), X0) = "(Ц+f[^+X™(3)]; a-(2)=x;«(2)+х™(2); х,(3) = X">(3)-f[Xi"(1)-
X,(4) = Xj"(0)-X»>(0); X,(5) = X<"(1)X<!>(1) + X«>(3)];
X(6) = X")(2)-X"(2); X,(7) = X™(3)-х(-1(3)].
Итерация 2 (m = 2, k = 0, 1):
^(О) = л\; (0) • А ЛО): Х«\1) = Х?\1) + Х™(1У, Х*\2) = .\7 (0) .V '.(О): ^>(3) = -V (1) (1 )-
где
■V ; (0) = д-(у') + л-(./ - 4); .V (1) = хЦ) - х(] - 4).
(2)т -
Рис. 3. Сигнальный граф полного скользящего
преобразования Хартли по основанию 2 с прореживанием по времени для N=8
Сигнальный граф этого алгоритма представлен на рис. 3. Алгоритм содержит 2
умножения на вещественные константы л/2 / 2 и 16 сложений, что соответствует оценкам (22) и (23) при #=8.
Сигнальные графы скользящих преобразований Хартли по форме отличаются от графов БПФ Кули-Тьюки и содержат дополнительные узлы. Зато все выполняемые в них операции вещественные.
Приведенные результаты, включающие свойства базисных функций и методы синтеза быстрых преобразований, составляют теоретическую основу спектрального анализа в базисе Хартли. Его применение особенно перспективно при решении спектральными методами различных задач обработки действительных сигналов.
Настоящая статья подготовлена по результатам НИР в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России на 2009-2013 годы (Государственный контракт № П1264 от 27.08.2009 г.).
СПИСОК ЛИТЕРАТУРЫ
1. Залманзон Л.А. Преобразования Фурье, Уолша, Хаара и их применение в управлении, связи и других областях. - М.: Наука, 1989. - 496 с.
2. Брейсуэл Р. Пребразование Хартли: Пер. с англ. - М.: Мир, 1990. - 175 с.
3. Трахтман А.М. Введение в обобщенную спектральную теорию сигналов. - М.: Сов. радио, 1972. - 345 с.
4. Трахтман А.М., Трахтман В.А. Основы теории дискретных сигналов на конечных интервалах. - М.: Сов. радио, 1975. - 208 с.
5. Власенко В.А., Лаппа Ю.М., Ярославский Л.П. Методы синтеза быстрых алгоритмов свертки и спектрального анализа сигналов. - М.: Наука, 1990. - 180 с.
6. Сюзев В.В. Быстрые преобразования Фурье для скользящего анализа спектра // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. - 1998. - N 2. - С. 29-38.
