Рекурсивний метод частотно-фазового аналізу двохтонального частотно-маніпульованого сигналу Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 621.372.542

І.А. БЕЗВЕРБНИЙ

РЕКУРСИВНИЙ МЕТОД ЧАСТОТНО-ФАЗОВОГО АНАЛІЗУ ДВОХТОНАЛЬНОГО ЧАСТОТНО-МАНІПУЛЬОВАНОГО СИГНАЛУ

Abstract: The two-tone frequency-shift keyed signal frequency and phase analysis technique by means of digital processing is proposed, the computation process using digital-analytical method on the base of slide spectral analysis is considered, the experimental check results of the tone signal program-demodulator system efficiency is described in this paper.

Key words: frequency and phase analysis, two-tone signal, numerical and analytical procedure, sliding spectral analysis.

Анотація: У роботі пропонується метод частотно-фазового аналізу двохтонального частотно-маніпульованого сигналу засобами цифрової обробки, розглядається обчислювальний процес за допомогою чисельно-аналітичного методу, побудованого на основі ковзного спектрального аналізу, описані результати експериментальної перевірки ефективності роботи програми-демодулятора тонального сигналу.

Ключові слова: частотно-фазововий аналіз, двохтональний сигнал, чисельно-аналітичний метод, ковзний спектральний аналіз.

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

Ключевые слова: частотно-фазовый анализ, двухтональный сигнал, численно-аналитический метод, скользящий спектральный анализ.

1. Вступ

Частотно-маніпульований сигнал складається з відрізків сум двох гармонічних коливань. Носієм інформації в сигналі є випадковий набір конкретних частот. Довжина відрізку може бути стандартною або визначатися за допомогою обертання фази. Найпоширенішим способом застосування цифрових частотно-маніпульованих двохтональних сигналів є телефонний зв'язок з тастатурним викликом у вигляді “двох з восьми можливих частот” - формат DTMF (Dual-Tone MultiFrequency). В сучасних системах розпізнавання такого сигналу відбувається за допомогою детектора, що складається з каскаду режекторних і смугових фільтрів. Докладний розгляд такого детектора міститься в [1, 2]. Поряд з тим двохтональний сигнал може бути використаний також в телеметричних, радіо- і гідроакустичних системах і системах космічного зв'язку.

Використання частотно-маніпульованого багатотонального зв'язку не обмежується форматом DTMF. В акустичних системах застосовуються аналізатори тональних сигналів, у яких передбачений суттєвий відхід частот сигналу. Тобто тони сигналу можуть мати майже довільну частоту, обмежену лише визначеним діапазоном. Наведений у статті наближений метод скороченого розпізнавання частот двохтонального сигналу може бути використаний в системах, побудованих за таким зразком. Алгоритм розпізнавання необхідних частот не потребує цифрової фільтрації, тому зменшує кількість і відповідно час обчислень. Крім того, він дозволяє рознести частоти в одному маніпуляційному відрізку в інших комбінаціях і в такий спосіб збільшити кількість сигналів у частотній смузі.

2. Сутність проблеми

Апаратурна фільтрація сигналу, що нині використовується все менше, висуває певні вимоги до габаритів АТС, тому на основі програмованого процесора цифрових сигналів розроблені

різноманітні цифрові варіанти тастатурних приймачів викликів. Такі системи є програмами, що складаються з багатосекційних фільтрів і передбачають велику кількість обчислень і великий об'єм пам'яті. Інтенсивне використання ресурсів цифрової апаратури в свою чергу призводить до високої ціни за обладнання, яке працює за даною схемою. Йдеться про роботу окремих апаратних модулів малого розміру, де ресурсна база буде завжди залишатися обмеженою. Пошуки нових методів частотно-фазового аналізу продиктовані метою зменшення кількості обчислень і у такий спосіб здешевлення апаратури.

3. Мета дослідження

Найважливішим завданням тастатурного аналізу двохтональних сигналів є частотна демодуляція тонів. Допоміжну інформацію може нести також початкова фаза кожного з маніпульованих тонів. Метою роботи є розробка методу швидкісного обчислення частот сигналу, про який відомо, що його спектр налічує дві частоти. Такий метод дозволить отримати дуже точні значення початкової фази і частоти обох маніпульованих тонів за допомогою декількох рекурсивних обчислень, реалізованих на дешевих сигнальних процесорах для роботи в реальному часі.

4. Ступінь розробки проблеми

В роботах [3, 4] докладно розглянуто чисельно-аналітичний метод, що передбачає визначення частоти і фази однотонального частотно-маніпульованого сигналу. Розглянемо принципові положення методу. Для визначення спектру за допомогою ДПФ необхідна кількість вибірок

маніпульованого тону рівна 2м +1. Числовий ряд спектра, згідно з означеним методом, дозволяє отримати точне значення початкової фази і частоти такого сигналу за допомогою нерекурсивної формули. В основі цього методу лежать принципи ковзного спектрального аналізу [5] і чисельно-аналітичний метод частотного аналізу однотональних сигналів [6]. У статті міститься розширення цього методу для двохтонального частотно-маніпульованого сигналу. Це передбачає побудову системи чотирьох рівнянь з чотирма невідомими. Її можна розв'язати за допомогою чисельного методу вирішення системи нелінійних рівнянь.

5. Фазовий аналіз

Вхідний двохтональний сигнал розглядається у вигляді

сигналу. Пошук значень /¡вх і /вх розглядається нижче. В даному випадку вважаємо їх відомими. Як відомо, фазовий спектр тонального сигналу можна розглядати як масив:

(1)

де /вх, /вх - частоти першого і другого тонів сигналу, в\, #2 - фази першого і другого тонів

(2)

де і - номер гармоніки; (М)) - дійсна частина гармоніки; (М )) - уявна частина.

Методика дослідження полягає в отриманні виразу (2) аналітичним шляхом. З цією метою на N

N -1 N -1

вибірок накладається прямокутне вікно від----------до -------- одиничної функції. Тоді спектральна

2 2

характеристика оов-перетворення, що відповідає першому тону сигналів типу /(п), подається

формулою

N-1 2

Яе(^ (М)) = Ке(Г (/І, /2,4/і, 42)) = 2 / (п)ссвС/1 п)

N -1

N-1 2

2 (с°8(/1ВХ п + & )+ С08(/2ВХ п + &2 )) С08(/1 п)

N-1

ґ N-1 2

2С08(/Г п)С08(/1п)

N -1 п=--------

V 2

С08

&)

ґ N-1 2

]Т8Іп(/1ВХ п) С08(/ п)

N-1 п=------

V 2

8Іп

(& ) +

ґ N-1 2

+

2С08/2ВХ п)С08(/1 п)

N-1 п=------

V 2

С08

ґ N-1 2

28Іп(/2ВХ п) С08^^/1 п)

N-1 п=------

V 2

8Іп

&2 ),

(3)

і відповідно спектральна характеристика віп-перетворення:

N-1 2

(М)) = Іт(7 (/1, у2,4/1, )) = 2 / (п)8Іп(/1 п) =

N-1

N-1 2

= 2 (с08(/1ВХ П + & )+ С08(/2ВХ П + &2 )) 8Іп(/1 п) =

N -1

ґ N-1 2

2С08/Г п)8Іп(/1п)

N-1

V ~2~

ґ N-1

С08

(& )

ґ N-1 2

28Іп/1ВХ п)8Іп^./1п)

+

2С08/2ВХ п)8Іп(/1 п)

N-1

V " ~2

N-1

V ~2~

ґ N-1

8Іп

(&1 ) +

С08

& )■

28Іп (/2ВХ п)8Іп(/1 п)

N-1

V " ~2

8Іп

&2 ), (4)

п=

2

п=

2

п=

2

п=

2

де /1 - частота функції, що належить до базису перетворення Фур'є і відповідає першому тону сигналу (1). Запроваджуємо наступні позначення:

N-i 2

N-1 2

A(fi )= Zcos(fiBX n )cosC/in) - B(fi )= Zsin(fiBX n )cos(fin) - (5- 6)

N-1 2

N-i 2

_ N-i ' 2

N-i 2

C (fi)= ZcosfiBX n) sin(fin) - D(fi)= Zsin (fiBX n) sin(fin) - (7- 8)

N-i 2

N-i 2

__ N-i ’ 2

N-i 2

^C/i)= ]Ccos(f2BX n) cos(fi n) - X(fi )= Zsinf2BX n)cos(fi n) - (9- 10)

N-i 2

N-i 2

N-i 2

N-i 2

7C/i )= ]Ccos(/2BX n) sin(fi n) - Z(fi )= Zsin(/2BX n) sin(fi n) - (11- 12)

N-i

N-i

С(#2 ) = Яе(^ (і® )) - и(/і ) ■ СОз(#2 ) - Х (Л ) • 5Іп(й2 ).

Н(#2 ) = Яе^(і® ))- Г(/і )• СОз(#2 )- г(/і ) ■ яп(б>2 ),

РІЙ / )= С(#2 ) = Яе(^ (І® ))-и (/і ) • сО^(#2 )- Х (/і )• ап(й2 )

' 2> М Н (Й2 ) Ьп^ (і® ))-Г (/і ) • ЄО8(Й2 )-г (/і )• 5Іп(Й2 ) '

де 5 - номер першого тону сигналу. Тоді вираз (2) може бути представлено

А{/і) • СОВ# )- В(/і )• 8Ш(йі )

С. (/і ) • СО8(йі ) - 0{/і ) • 8Ш (йі ) .

А(/і)-С (/і )• РЙ2, /і) вЦі)-в(/і )• Р(Й2, /і) ■

Таким чином, отримано значення # через #2 . Аналогічним шляхом будується ще одне рівняння, за допомогою якого буде отримано значення #2 через #. В результаті аналогічних елементарних перетворень нове рівняння буде повністю повторювати попереднє, але значення індексів 1 і 2 зміниться 2 і 1 відповідно. Воно має такий розв'язок:

А(/2)-С (/2 )• Р#і, /2)

Р(вг, fi ) =

Рівняння має такий розв'язок:

q = arctan

(13)

(14)

(15)

(16)

(17)

b(/2 )-d(/2 )• p(q, /2)’

(18)

де

n=

n—

n=

n—

n—

n=

2

2

N-1 N -1

2 t \ 2

г-вх „ І_______í r \ r>í j- \ V _• I .свх

A(f2 )= n)C0s(f2 n), B{f2 )= ]^SÍn(/2BX П) C0s(f2 n) , (19, 20)

N -1 N-1

n=- n=-

22

N-1 N-1

2 1 г-вх Л ■„( r \ n// ^ • I /вх

C(/2 )= ]ССте(/2вх n) sin(f2 n) , D(/2 )= ZSin(/2вх n) SÍn(/2 n), (21, 22)

N-1 N-1

n=- n=-

22

p(0 f )= Re(F(ww ))- U(/2 )- C0s(02 )- X(/2 )- Sin(q2 ) (23)

0 fl) = Im(F(jws))-Y(/2 )■ C0S&)-z(/2 )■ sin(q) , ( )

де s - номер другого тону сигналу,

N-1 N-1

2 вх 2 вх

U(/2 )= n) C0s(f2 n), X(f2 )= ZSin(/1вх n) C0s(f2 n) , (24, 25)

N-1 N-1

n=-- n=--

22

N-1 N-1

Y (/2 )= n) Sin (/2 n) , Z (f2 )= ZSin (/1вх n) Sin(f2 n) ■ (26, 27)

N-1 N-1

n=-- n=--

22

Маємо систему з двома невідомими 01 і 02 , що складається з двох рівнянь. Вважаємо /1вх

і /вх в такому випадку відомими. Насправді їх визначення потребує ще двох рівнянь. Побудова

цих рівнянь розглядається нижче.

6. Частотний аналіз

Частотний спектр не залежить від фазового, тому для спрощення розглядається сигнал з нульовою фазою в кожному з тонів. У такому випадку вхідний сигнал має вигляд функції f (n):

f (n)= C0S(/[вх n)+ C0s(f2х n)= C0S2p(f1 +A/1)n + C0S2p(f2 +Af2)n , (28)

де /г, /2 - ціла складова частоти першого і другого тонів сигналу, А/г, А/2 - дробова складова

частоти першого і другого тонів сигналу, п - номер вибірки, ціле число, 0 < п < N , N - непарне число, кількість вибірок на період аналізу сигналу.

Подібно до методики визначення фази, у випадку визначення частоти необхідно отримати аналітичний вираз для обчислення спектральної характеристики, що відповідає першому тону. В результаті накладання прямокутного вікна отримуємо аналітичний вираз спектральної характеристики перетворення Фур'є для першого тону функцій вигляду /(п) :

F CM) = Y (fi, f2, Afl, A/2) =

N-1 2 -j

2k/1n I f (n)e 'y N

N-1

N-1 2

I

__ N-1 _ 2

cos

2k / n 2k 4/ n

f

N

- +

N

+ cos

2k f2 n 2kA/2 n

V

N

- +

N

- / 2Kf1> , \ N

n_

2

V

sin(p4f1 )

2'SinfPr'

+ ■

sin (кА/ )

_ . f2-k-f1 K-Afl 2-sinl-----^ + J1

■ +

N

N

sin

(k4/2 )

sin

(k4/2 )

+----------7-------------------------'T +-------/

_ . fpf2 kf1 k4/2^ . fpf2 kf1 k4/2

2 - sinl -----------------------------------------------— + —— I 2 - sinl + - 2

NN

N

NN

N

Визначимо

G( m)_ g(/1, f2, a/ )_ “KL +-------------------------------------2^

h ( M )_ h ( f,, f2, a/2 ):

. k4/1 . f 2k f1 kA/ ^ ’

2- sin—^ 2-sinl—I

N ^ N N )

sin(p4f2 )

sin(p4f2 )

+ ■

2 sin

+--+ -

NN

Тоді |F (jM )) _ |G (f,, f2, 4/i ) + H (f,, f,, 4f2 ).

З іншого боку,

|F ( j(w, +1)) _ Y [( f, +1), (4/1 -1), 4/i, 4/2 11

N-1

”2”

I

N-1

cosf 2p(f1 +!) n + 2k(4/1 - !) n ^ +

N

N

2

+ cos

f 2k(/2 +1) n + 2k (4/2 -1)

V

N

+ -

n

N

/ 2к(/і+1) n

N

N

(29)

(30)

(31)

(32)

n_

e

sin(p 4/1 )

2-sinfKM

l N

+ ■

sin(p 4/1 )

K

N.

_ . f 2k f kAf k

2 - sinl-----^ ^ +

- +

N

N N.

Йй (ж А/ )

.. (ж /2 ж / жА/2 ж

2 - БІПІ ——---------------------— + - 2-

БІЙ

(жА/2 )

N N

N N.

„ . (ж /2 ж / жА/2 ж

2 - БІйІ + —— + —

І N N N N.

Визначимо

0(] (ю,+1)) = 0( /, +1, /2, А/, ) =

Бій (ж А/1)

- + -

БІЙ (ж А/,)

(34)

V N N) V N

Н (](щ + 1)) = Н (/, + 1, /2, А/2 ):

БІЙ (жА/2 ) БІЙ

N N.

+ -

(жА/2 )

_ . (ж /2 ж / жА/2 ж ї _ . (ж /2 ж / ж А/2 ж

2 - БІЙІ ---— + ——-----І 2 - БІЙІ +-— + —

V N NN N) V N N NN у

Тоді ^ (](щ +1)) = \о (/, +1, /2, А/і)+Н (/1 +1, /2, А/2).

(35)

(36)

Нехай Л(ю і) і Л(ю і+1) - відомі значення гармонік, отримані за допомогою ДПФ. Розглянемо наступне відношення:

К(т)= К(/вх)= Л((°і )-Н(/1,/2,А/2) = ^(/1,/2,А/1) (37)

К^2 )= Л(рм)-н(/1 +1,/2,а/2) = о/ +1,/2,АЛ)• (37)

Останнє співвідношення повністю співпадає зі співвідношенням, отриманим для одного тону, що докладно розглядалося у [4]. Відповідно вирішення такого рівняння не відрізняється від отриманого раніше, і вхідна частота у такому випадку знаходиться у вигляді

N Л(/“)

/1 = — агйЯЙ

ж

В

(38)

де

л(/“ )= -К (/“) -

БІЙ2 ж/1 + К(/,х )-

N

• 2 ж - /і . 2 ж - (/і + 1)

БІЙ2-------— - БІЙ2---------—--------- +

2 ' N N

ж - /1 • 2 ж - (/1 + 1) ж - (/1 + 1)

+ СОБ-----БІЙ2----—-------СОБ-

N

N

N

(39)

В

(/вх ) = -К(/вх ) - БІЙ2 + К(/вх ) - БІГ 2 ж - /1 — 2 ж - (/1 + 1)

N

БІЙ

N

+

ж- /1 -2 ж - (/1 + 1) ж - (/1 + 1)

+ СОБ-------—- БІЙ2------------------------------------------—- - СОБ-—- +

N

N

N

+К

(/“ )к(/“ )- БІЙ2

ж-О! +1) _ж-/1 л-(/1 +1)

СОБ

СОБ -

(40)

N N N

Таким чином, отримано значення /вх через /вх. Аналогічним шляхом будується ще одне

„вх _р вх

рівняння, за допомогою якого буде отримано значення / через /1 . В результаті аналогічних

перетворень нове рівняння буде повністю повторювати попереднє, але значення індексів 1 і 2 зміниться 2 і 1 відповідно. Вхідна частота в такому випадку знаходиться у вигляді

f2ВХ = ~ -arctan p

"Де]

f)'

(41)

де

A(f вх )= -K(fвх) - sin2 pf2 + K(f вх)- sin2 pf2- sin2 p'f 1 +

VI ! \J I ! N v 1 ! N N

p- f 2 • 2 p - (f 2 +1) p- ( f 2 +1)

+ cos—--sin2—^------- - cos——------------------------- , (42)

N N N

B(f вх )= -K(fВХ) - sin2 + K(f ВХ)- sin2 - sin2 p'f 1 +

V-7 \ ) V x ) N v 1 ! N N

p - f2 • 2 p - (f2 + 1) p - (f2 + 1)

+ cos —— - sin2 ——--------------------- - cos-—-- +

N N N

+ K (f ВХ )-K (f ВХ) - sin2 p-(f2 +1) - cos - cos p-(f2 +1) , (43)

1 1 N N N

де K (fBX) - відношення гармоніки другого тону і сусідньої.

Отже, маємо систему, що складається з чотирьох рівнянь. Систему найлегше вирішити шляхом послідовних наближень за допомогою „процесу Зейделя” [7]. Цей метод має такий символічний вигляд:

I *„+1 =j(xn, Уп )

I Уп+1 = j (*п+Ь Уп ) .

Точність методу досить висока і в нашому випадку дає результат з машинною похибкою після максимум 10 ітерації. Розв'язок методом Зейделя наведених вище систем є в багатьох відношеннях кращим за розв'язок за допомогою метода Ньютона, бо останній вимагає визначення похідних, що може ускладнити вираз і у такий спосіб бути додатковою перешкодою у випадку необхідності запрограмувати метод мовами низького рівня.

Поданий метод аналізу двохтонального сигналу має певні обмеження, які відрізняються від обмежень, що виникають підчас аналізу однотонального сигналу. Завдяки тому, що метод є наближеним, А/ та А/2 може дорівнювати асимптотичному значенню 0,5. Крім того, обидва тони

можуть мати одну і ту саму грубо визначену частоту. Наприклад, метод дозволяє визначати сигнал, що містить дві частоти 2,6 Гц і 2,8 Гц. Натомість виникають обмеження наступного характеру. Грубо визначені частоти обох тонів повинні бути або обидві парні, або обидві непарні. Інакше інтерполяційний ряд буде розходитися.

7. Практичне застосування методу

Нехай на вхід демодулятора надходить сигнал, що складається з двох тонів:

2 p-2-87Гц »0282рад;

64 ’ сек

2-p-6,62 Гц » 0649 рад.

64 сек

Відбувається квантування сигналу з частотою 6 кГц, тобто 65 вибірок надходять за період

0,0107 секунди. Вигляд відцифрованого сигналу зображено на рис. 1. З масиву з 65 вибірок знаходиться масив ДПФ (рис. 2), що в результаті обробки дозволяє отримати результати, представлені в табл. 1. В даному випадку точне значення може бути отримане за 4 ітерації.

Таким чином, для передачі маніпульованого сигналу достатньо 65 вибірок. Обертання фази на 90° може розглядатися як сигнал синхронізації першого маніпульованого відрізка. Наступні маніпульовані відрізки синхронізації не потребують. В такий спосіб за 1 секунду може бути передано близько 100 двійкових комбінацій частот. Кількість тонів у таких комбінаціях, як можна бачити з можливостей методу, дозволяється не обмежувати числом 8. Метод передбачає, щоб сигнал складався лише з двох тонів, але він може бути розширений до роботи з більшою кількістю тонів. При цьому збільшення частоти дискретизації не обов'язкове.

8. Висновки

Наведений метод випробувано програмно на сучасному комп'ютері, оснащеному процесором класу Pentium!I, на операційній системі Windows98 за допомогою програмного забезпечення Cool Edit 2000, Borland Delphi версії 6.0 фірми Enterprise і бібліотеки WinAPI і дозволяє з дуже високою точністю визначати частоту і початкову фазу вхідної функції частотно-маніпульованого сигналу.

Безперечною перевагою даного методу є мінімізація обчислень порівняно з вищезгаданими методами обробки М-сигналів, унаслідок чого створюється можливість обробки в реальному часі цифрових тональних сигналів з дуже близьким рознесенням частот. Слід зважати, що в результаті реалізації методу на процесорах, що не мають системи операцій з плаваючою комою, розмір програмного забезпечення може зростати, і разом з ним зростатиме час на роботу програми обробки сигналу. Особливо це стосується обчислення тригонометричних функцій. Таким чином, система частотно-фазового аналізу, розроблена на базі процесора ADSP21xx фірми Analog Device, дозволяє отримувати дані про зміну частот і початкових фаз вхідного сигналу за приблизно 8 мксек.

Наприклад, сигнальний процесор ADSP2188, максимальна частота якого 75 МГц, виконує одну операцію за 13 нс. Обчислення швидкого cos-перетворення займає близько 100 операцій.

0 5 10 15 20 25 30 35 40 45 50 55 60 Рис. 1. Вхідний сигнал в часовій області

1

1 1

1 1 1 1

.1 і ; і і І,

І J Шк, і4і] І II

0 20 40 60

Рис. 2. Вхідний сигнал в частотній області

Таблиця 1. Послідовна робота ітераційного метода

№ ітерації Перший тон Другий тон

1 2,7148 6,6389

2 2,8765 6,6190

3 2,8697 6,6201

4 2,8700 6,6200

Обчислення апроксимації косинуса, синуса і арктангенса займає відповідно по 20 операцій. Визначення косинуса і синуса під час обчислення частоти необхідно провадити 15 разів, згідно з формулами (38) і (41). Таким чином, визначення частоти без швидкого сов-перетворення виконуватиметься за близько 450 операцій. Обчислення фази реалізується за близько 250 операцій.

Для реалізації даного методу на практиці бажано, щоб сигнал не був спотворений адитивним шумом. Розгляд роботи методу з тональним сигналом і адитивним шумом виходить за межі теми статті. Але робота з таким сигналом реальна, хоча вимагає підключення підпрограми, яка містить близько 300 операцій на сигнальному процесорі ДР8Р2188.

СПИСОК ЛІТЕРАТУРИ

1. Семенов Н. MATLAB для DSP. Моделирование приемника DTMF. - Часть 1. -

http://chipnews.gaw.ru/html.cgi/arhiv_i/00_02/.

2. Семенов Н. MATLAB для DSP. Моделирование приемника DTMF. - Часть 2. -

http://chipnews.gaw.ru/html.cgi/arhiv_i/00_05/.

3. Безвербний І.А. Удосконалення чисельно-аналітичного методу спектрального аналізу тональних сигналів:

Зб. наук. пр. «Засоби комп'ютерної техніки з віртуальними функціями і нові інформаційні технології». - К.: Ін-т

кібернетики, 2002. - Т. 1. - С. 105-109.

4. Безвербний І.А. Чисельно-аналітичний метод демодуляції тональних сигналів з урахуванням прямокутного вікна // УСіМ. - 2005. - № 4. - C. 19-29.

5. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. - М.: Мир, 1978. - 848 с.

6. Семотюк М.В. Численно-аналитический метод спектрального анализа тональных сигналов // УСіМ. - 2001. -№1. - C. 36-42.

7. Демидович Б.П., Марон И.А. Основы вычислительной математики. - М.: Наука, 1966. - 664 с.