Быстрое дробное преобразование Фурье

Степан Александрович Мартюгин; Сергей Владимирович Поршнев

С. А. Мартюгин, С. В. Поршнев

Аннотация—Дробные преобразования Фурье (ДрПФ) являются семейством унитарных преобразований

/р "\2л . При а = л/2 ДрПФ принимает облик

I )а=0

классического преобразования Фурье ( р *=л/2 = Р ), при а = л становится инверсным преобразованием (р а=* = J ), при а = 3л /2 соответствует обратному

преобразованию Фурье ( F

= F 1 ), а при а = 0 (или

а = 2л) вырождается в тождественное преобразование (Fa=0 = I )• Семейство (f формирует

I J1X1

однопараметрическую непрерывную унитарную группу с аддитивным умножением F " • F р = F а+р. В аналитическом представлении ДрПФ является разложением сигнала по базису, состоящему из набора сигналов с быстро изменяющейся частотой (ЛЧМ-сигнал, или chirp).

Несмотря на то, что на сегодняшний день ДрПФ нашли широкое применение в различных задачах обработки сигналов и изображений, проблема совершенствования их быстрых алгоритмов остаётся актуальной. Большинство из представленных алгоритмов не могут обеспечить одновременно высокую точность вычислений и высокое быстродействие (сопоставимое с быстрым преобразованием Фурье (БПФ)). Целью работы является разработка устойчивых быстрых алгоритмов (БА) дискретных дробных и четырёх-параметрических преобразований Фурье (ДЧППФ) с минимальной вычислительной ошибкой. Для этого с помощью метода проективной декомпозиции оператора дискретного преобразования Фурье (ДПФ) были впервые получены аналитические выражения матричных элементов, составляющих ядро этих преобразований. Проведенный анализ вычислительной сложности полученных алгоритмов показывает, что их сложность сопоставима со сложностью быстрого преобразования Хартли (БПХ), O(N log N) ).

Ключевые слова—Обработка сигналов, дискретное преобразование Фурье, дробное преобразование Фурье.

I. Введение

Идея о дробных степенях оператора Фурье появилась в математической литературе в работе Е. Кондона 1937г. [1]. Позже, в 1961 году В. Баргманн в [2] дал более точное определение этого преобразования, основанное на многочленах Эрмита.

Статья получена 24.10.2023.

Степан Александрович Мартюгин, Уральский федеральный университет имени первого Президента России Б. Н. Ельцина, stmart2608@gmail.com

Сергей Владимирович Поршнев, д.т.н., профессор, Уральский федеральный университет имени первого Президента России Б. Н. Ельцина, ведущий научный сотрудник, Институт математики и механики им. Н.Н. Красовского Уральского отделения РАН, s.v.porshnev@urfu.ru

В 1980г Намиас заново открыл ДрПФ [3]. В своей работе, посвященной использованию дробного ПФ для решения некоторых задач, связанных с квантовым гармоническим осциллятором, он представил его операционное исчисление. Серьёзное развитие это направление получило в 90-е, после ряда публикаций Мендловича и Озактаса [4]-[6], которые предложили использовать ДрПФ в оптике. Значительный вклад в разработку теории дробного преобразования Фурье внесли Дикинсон и Штиглиц [7], Пей [8]-[9], Бултхил и Мартинез [10], Лабунец [11]-[12].

Сегодня ДрПФ находят широкое применения для решения различных задач в области обработки сигналов и изображений: фильтрации [13]-[14], сжатия [15], распознавания образов [16], квантовой механики [17], нейронных сетей [18]-[20] и др. [21].

Одной из фундаментальных проблем дискретных ДрПФ (ДДрПФ) является разработка быстрых алгоритмов их реализации. Озактас в своей работе [22] предложил способ вычисления ДрПФ, основанный на свёртке сигнала с ЛЧМ-функцией (чирп-функцией) ядра ДрПФ, полученной Баргманном в [2], [10], [12]. Вычислительная сложность этого преобразования совпадает со сложностью быстрого преобразования Фурье (БПФ), однако полученное данным способом преобразование не является коммутативным и инверсным, поэтому оно не подходит для решения практических задач [10].

Кандан [23], Дикинсон и Штиглиц [7], Грунбаум [24] предложили вычислять ДрПФ методом подбора коммутирующей с оператором ДПФ трёх-диагональной матрицы [25]. Известно [10], [25], что если оператор F коммутирует с оператором H (т.е. выполняется FH = HF), то они имеют общий набор собственных векторов. Таким образом, задача сводится к построению собственных векторов оператора H . Этот метод позволяет вычислять спектр ДрПФ с высокой точностью, но имеет большую вычислительную сложность O(N2) [10].

В статье проводится построение ядра ДДрПФ и нового дискретного четырёх-параметрического преобразования Фурье (ДЧППФ) с помощью проективной декомпозиции оператора ДПФ, после чего разрабатываются БА их точного вычисления. Показано, что вычислительная сложность разработанных алгоритмов совпадает со сложностью БПХ (O (N log N),

совпадает с БПФ).

Работа построена следующим образом: в секции II приведён краткий обзор теории ДрПФ, в секции III выполняется проективная декомпозиция оператора

ДПФ, секции IV и V посвящены построению ядра и БА ДДрПФ и ДЧППФ соответственно.

II. Дробное преобразование Фурье

Известно [2], что собственными функциями оператора преобразования Фурье являются функции Эрмита: 1

III. Проективная декомпозиция оператора ДПФ

Y ( х ) =

42"п\4ж

H„(. х)г-

jn

где И„ (х) = (_1)п e"1 j^e*

(1)

- многочлен Эрмита

порядка п е N, и -t

F [Y- (х)] = Т J Yя (х) e-^dx

(2)

= VY, (У) = e J 2П Yn (y),

где Хп = (-])" = е 2 является собственным числом,

соответствующим п-ой собственной функции. Тогда оператор ДПФ имеет следующее собственное разложение:

Г ( * м

F = U sdiag

■и-1

(3)

= £ e J 2 Yn (x) Yn (y),

где и = (/) | ^ (/) - матрица, составленная из

собственных векторов оператора ДПФ. Определение 1. Дробное преобразование Фурье получается возведением собственных значений оператора ДПФ в степень аеС :

Г ( * \\

F a = U ! diag

>и-1

да Л

- j—na

2

(4)

= Х Y n ( x ) Y n ( y ).

Существует также и а -параметризация, где а = ая 12, а ей . Оператор ДрПФ обладает следующими свойствами:

1) Периодичность с периодом 4 (2* в случае а -параметризации);

2) Коммутативность Б Т ь = Б ьБ ";

3) Инверсия Б " Б = и Л"-" и-1 = Б 0 = I ;

4) Сохранение энергии \/(х)|2 = " {/ (х)}| ;

5) Согласованность с тождественным оператором

р "=о = р а=2* = I и оператором ДПФ

р "=1 _ р а=*2 _ р .

В настоящее время дискретные аналоги функций Эрмита (которые являются собственными векторами оператора ДПФ) неизвестны [10], [23], [25], поэтому для построения ядра ДДрПФ и его быстрого алгоритма будет использована техника проективных операторов.

Пусть F =

2л - J—nm

с N

- оператор ДПФ, размером

N х N . Согласно спектральной теореме для унитарных операторов [26] это преобразование можно представить в виде суммы проекторов:

F = P

(5)

где одномерные действительными

проекторы P; симметричными

являются

( P = pa

идемпотентными (р2 = Р ), ортогональными (рр = 0 при к ф I), инвариантными (рр = \р) операторами.

Поскольку Р4 = I, то оператор Р имеет четыре различных собственных значения, которые являются

4 1

корнями характеристического уравнения 2 = 1, т.е.

* г

-к I \к

4 = е 2 =(-]) , к = 0,1,2,3. Следовательно, выражение (5) можно представить следующим образом:

3

р=Хл р

к =0

= 4 Р + 4 Р + 4 Р + 4 Р (6)

=р - ]р - р+;р.

Найдём выражения проекторов. Для этого подействуем Р на правые и левые части выражения (6) и запишем результат в матричной форме:

(7)

F0 "1111] Po

F1 1 -J -1 J P1

F2 1 -1 1 -1 P2

F3 ; J -1 -J. P3

Откуда

P = 4 ( I + / F + (/ F )2 +(j

(8)

(9)

где I - единичная матрица. Таким образом Ро,2 = 1 (I + J ± 2С) ,

Р1,3 = 1 (1 - J ± 28 ) ,

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

где С и 8 - матрицы дискретного косинусного и синусного преобразований (ДКП и ДСП соответственно) порядка N, элементы которых

[с] =-1

1 ^ у/Ы ,

( 2л cos — mn

l N

[S] = "1

L Jnm JN

' . ( 2л sin — mn

l N

1 0

0 0

0 1 •■• 0

IV. Быстрое дробное преобразование Фурье

Определение 2. ДДрПФ является взвешенной суммой четырёх проекторов:

N-1

k=0

71

n=0

N-1

n n n

- j—a - j—2a - j—3a

2 2 2

¥" = р + е 2 р + е 2 р + е 2 Р3 (10)

Найдем выражение ядра этого преобразования, выполнив замену а = ап /2 и используя (9):

Г = р + е~ ар + е~3 2ар + е~3 3аР3

= 1 е"30а (I + J + 2С) +1 е"а (I - J + 28)

+1 е"32а (I + J - 2С) +1 е"33а (I - J - 28 ) (11)

а

XN =

V2

IM-i V2 -IM-i

IM-i I M-i

V2

Im -Im

Im IM

нечетного N , I =

i .

матрица.

Прямые вычисления показывают, что

xn fn xn =

-,[ X]

-JSlX]

XNCNXN i XN (FN + FN ) XN

XNSNXN = - J i XN (fn - FN ) XN =

'/-[ X]

CM+i

G

S[X]

T V __Ü2 V _

XN J N XN XN FN XN

Подставляя (13)—(15) в (12), получаем

YN TTа Y

XN FN X N

(

+ e j sin«

ja/2T

e i

-jЗа/2

a а

cos —1„ + i sin — 2 N 2

-i[ X]

1M+i л

IL _ У

GM+i л

S[X]

_ У

- J—/2x

e J I,

= e iae ia/2 cos a I cos — I + i sin — J 1 2 2

+ e" iai sin—(C + e"ia S ).

Для построения быстрого алгоритма ДДрПФ воспользуемся свойством симметрии матриц ДПФ. Известно, что матрицы ДПФ без первого столбца и первой строки являются центрально-симметричными. Поэтому они приводятся к блочно-диагональному виду следующими матрицами соответственно для четного и нечётного N :

V2

(

cos—I^ + J sin»

Отсюда найдём Fa, принимая во внимание, что

-I

C[X] !

; s[x]

X,

X]

sLx]

XN = Нд, - преобразование Хартли:

F" = e~3j"/2 |cos "i + jsin cJcos"I + jsincH). (16)

Оценим вычислительную алгоритма. Пусть s -вектор сигнальных отсчётов. Тогда его ДрПФ-спектром (5") будет:

S" =F"•s

а

Xn = T

где M = N /2 для четного N, M = ( N -1) /2 для

. i"

_ -3 Ja/2

= e-3j—/2

cos — I + J sin — J [•(cosal + J sinaH)• s cos — I + J sin — J |-(cosaI • s + J sinaH • s)

- антидиагональная

(12)

где - С®! и 8® - матрицы ДКП и ДСП соответственно, Ь = М — 1 для чётных N и Ь = М для нечётных N. Поэтому

(13)

(14)

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(15)

Следовательно, вычисление спектра Sa представляет собой последовательность следующих операций:

1) Вычисление суммы единичных матриц cos (a/2) Iw + i sin (a/2) J = X^ - сложность O (N).

2) Вычисление произведений cos a • IN • = in -сложность O (N).

3) Вычисление быстрого БПХ i sina^ • ^ = h -

сложность O (N log N) . Отметим, что пп. 1-3

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

4) Вычисление суммы векторов /я + кп = h -

сложность O (N).

5) Вычисление произведения X'Nh'n - сложность

O (N).

Таким образом, вычислительная сложность алгоритма ДДрПФ совпадает со сложностью наиболее вычислительно сложной операции - БПХ [27]: O (N) + O (N) + O (N log N) + O (N) + O (N ) =

= O (N + N + N log N + N+N ) = O (N log N).

V. Быстрое четырёх-параметрическое

ПРЕОБРАЗОВАНИЕ ФУРЬЕ

Определение 3. ДЧППФ будем называть следующее преобразование:

Л . Л Л г, Л _

Sr, г, г, г, \ - 4a0 - j—a - j — 2a7 - j—3 a

рН.^л} = eJ2 0p + eJ21P + eJ2 2p + e 2 3p (17)

2

Найдём ядро этого преобразования, используя (9):

Fio0.a1.a2.03} = e~J4ao р + e~Ja1 P + e j2"2 P + eJ3ai P

= 1 e~J400 (I + J + 2C) +1 e~Ja1 (I - J + 2S)

1

+ -e

J 2a2

1

(I + J - 2C) + - e-J3° (I - J - 2S )

функций ДДрПФ и ДЧППФ, на основе которых разрабатываются их быстрые алгоритмы. Полученные дискретные преобразования являются коммутативными и инверсными, а их БА имеют сложность O (N log N) .

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

= 1 (e~J400 + e Ja )(I + J) +1 (e~J4a° - e~J2a ) C +1 (e" J1a1 + e"J3a ) (I - J) + 1(e-J1a - e"Ja )S.

Обозначим ß = e -J 400 + e ~J 202

ß = e-J4a° -e-J2ai

(18)

г = е-п<ч и У = еч1а^ -еЧЪщ . Тогда:

Для построения быстрого алгоритма ДЧППФ построим X-представление матрицы р^л^л! (11)-(14):

XT F'"

XNF N

,0,a3}

X = 1

Xn = 2

ß Л

Л

П L

+ -

ß\

LM+1 _!_

Vp[ X]

CM+1

\nL

¡¡[ X]

Следовательно:

F'

l{a0.a1.a2.aз} _

= 1 ((ß + r)I + (ß-r) J )

+ -A((ß + f)\ + (ß-f)i)-\i (17)

Л(х'+х'-н),

где X' = (ß + y)l + (ß-y)3, X' =(ß + f)l + (ß-f)3. Вычисление ДЧППФ-спектра сигнала s :

S

= F

= I(X'+X'.H)., = -s + X' -H-s)

является последовательностью следующих операций:

1) Вычисление произведения

XN • sn = ((Р + у) In + (Р~у) Jn )• sn - сложность

2) Вычисление произведения X'N • Н • \д . -сложность

O (N log N) + O (N) = O (N log N).

3) Вычисление суммы полученных в п.1 и п.2 векторов - сложность O (N).

Таким образом, вычисление спектра быстрого ДЧППФ, как и в случае быстрого ДДрПФ совпадает со сложностью БПХ:

O (N ) + O ( N log N ) + O ( N ) = O (N log N ).

VI. Заключение В статье получены аналитические выражения ядерных

Библиография

[1] E. U. Condon "Immersion of the Fourier transform in a continuous group of functional transforms", Proc. Nat. Acad. Sci, 1937, Vol., 12, pp. 158-164.

[2] V. Bargmann "On a Hilbert space of analytic functions and an associated integral transform. Part 1", Commun. Pure Appl. Math., 1961, 14, pp.187-214.

[3] V. Namias "The fractional order Fourier transform and its application to quatum mechanics", J. Inst. Math. Appl., 1980, Vol. 25, pp. 241-265.

[4] H. M. Ozaktas. and D. Mendlovic "Fourier transform of fractional order and their optical interpretation", Opt. Commun., 1993, 110, pp. 163169.

[5] H. M. Ozaktas, M. A. Kutay, and D. Mendlovic, "Introduction to the fractional Fourier transform and its applications", Advances in Imaging Electronics and Physics. New York: Academic, 1999, ch. 4.

[6] H. M. Ozaktas, B. Barshan, D. Mendlovic, and L. Onural, "Convolution, filtering, and multiplexing in fractional Fourier domains and their relation to chirp and wavelet transforms," J. Opt. Soc. Amer. A., vol. 11, pp. 547-559, 1994.

[7] B. W. Dickinson and K. Steiglitz. "Eigenvectors and functions of the discrete Fourier transform" IEEE Transactions on Acoustics, Speech and Signal Processing, vol 30(1), pp. 25-31, 1982.

[8] S. C. Pei and M. H. Yeh, "Improved discrete fractional fourier transform," Opt. Lett., vol. 22, pp. 1047-1049, 1997.

[9] S. C. Pei, C. C. Tseng, M. H. Yeh, and J. J. Shyu, "Discrete fractional Hartley and Fourier transforms," IEEE Trans. Circuits Syst. II, vol. 45, pp. 665-675, 1998.

[10] A. Bultheel, and H. Martinez, "Computation of the Fractional Fourier Transform, preprint".

[11] E. Labunets, V. Labunets, "Fast fractional Fourier transforms", Proc. of Eusipco-98, Rhodes, Greece, 8-11 Sept. 1998, pp. 1757-1760.

[12] E. Ostheimer, V. Labunets, and S. Martyugin "Fast Infinitesimal Fourier Transform for Signal and Image Processing via Multiparametric and Fractional Fourier Transforms", Supplementary Proceedings of the 4th International Conference on Analysis of Images, Social Networks and Texts (AIST'2015), pp. 19-27.

[13] L. Chen, D. Zhao "Application of fractional Fourier transform on spatial filtering", Optik 117, pp. 107-110, 2006

[14] M. F. Erden, M. A. Kutay, H. M. Ozaktas, "Applications of the fractional Fourier transform to filtering, estimation and restoration" Proceedings of the IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing (NSIP'99), Antalya, Turkey, June 20-23, 1999, pp. 481-485.

[15] C. Vijaya, J. S. Bhat, "Signal compression using discrete fractional Fourier transform and set partitioning in hierarchical tree", Signal Processing 86 (8) pp. 1976-1983, 2006.

[16] A. W. Lohmann, Z. Zalevsky, D. Mendlovic, "Synthesis of pattern recognition filters for fractional Fourier processing", Optics Communications 128 (4-6) pp. 199-204, 1996.

[17] F. Hong-Yi, F. Yue, "Fractional Fourier transfromation for quantum mechanical wave functions studied by virtue IWOP technique" Commun. Theor. Phys. (Beijing, China) 39 pp 417-420, 2003.

[18] B. Barshan, B. Ayrulu, "Fractional Fourier transform pre-processing for neural networks and its application to object recognition", Neural Networks 15 (1) (2002) 131-140.

[19] X. WU, R. TAO, D. HONG, Y. WANG, "The FrFT convolutional face: toward robust face recognition using the fractional Fourier transform and convolutional neural networks", SCIENCE CHINA Information Sciences, Volume 63, Issue 1: 119103 (2020)

[20] F. §ahinuj, A. Koj, "Fractional Fourier Transform Meets Transformer Encoder", IEEE Signal Processing Letters, Vol 29, pp. 2258-2262, 2022

[21] S. Ervin, D. Igor, S. Ljubisa, "Fractional Fourier Transform as a Signal Processing Tool: An Overview of Recent Developments", Signal Processing 91 (2011) pp. 1351 - 1369.

[22] H. M. Ozaktas, O. Arikan, M. A. Kutay, and G. Bozdagi, "Digital computation of the fractional Fourier transform," IEEE Trans. Signal Processing, vol. 44, pp. 2141-2150, Sept. 1996.

1

2

[23] C. Candan "The Discrete Fractional Fourier Transform", M.S. thesis, Bilkent Univ., Ankara, Turkey, (1998).

[24] F. A. Grunbaum, "The eigenvectors of the discrete Fourier transform: A version of the Hermite functions", J. Math. Anal. Appl. 88(2) (1982) pp. 355-363.

[25] A. Serbes, L. Durak-Ata, "The discrete fractional Fourier transform based on the DFT matrix", Signal Processing 91 (2011) 571-581.

[26] М. С. Беспалов, "Спектральное разложение дискретного оператора Фурье", доступно по адресу: http://dha.spb.ru/PDF/DFTSpectrum.pdf

[27] Д. Кнут "Искусство программирования, том 1. Основные алгоритмы" 3-е издание, стр. 139.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Fast fractional Fourier transform

S. A. Martiugin, S. V. Porshnev

Abstract—Fractional Fourier transform (FrFT) is one-

parametric family of unitary transformations (f a\2n .

I Ja=0

FrFT are usually interpreted as rotations in the time-frequency plane. In case a = n /2 FrFT becomes an ordinary Fourier transform (F a=n/2 = F ), in case a = n we have inversion transform (F a=n = J ), for a = 3n /2 FrFT

corresponds to inverse Fourier Transform ( F

= F -1 ),

and for a = 0 (or a = ) we have identity transform (f 01=0 = i ). The family (f ") forms a one-parameter

I t aeD

continuous unitary group with additive multiplication F a • F p = F a+p. FrFT corresponds to decomposition of a signal into linear combination of chirps (or signals in which the frequency increases (up-chirp) or decreases (down-chirp) with time).

FrFTs have found wide application in various signal and image processing tasks, however the problem of improving their fast algorithms remains relevant. In this work we develop fast algorithms for discrete fractional and four-parametric Fourier transforms (FPFT) with minimal computational error. For that, we derive kernels for these transforms using a method of projective decomposition of the discrete Fourier transform operator (DFT). Analysis of computational complexity shows that the complexity of obtained algorithms is comparable to the fast Hartley transform (FHT, O ( N log N) ).

Keywords—Signal processing, Fourier Transform, Fractional Fourier transform

References