Розробка

Реалізація цілочисельного ШПФ на ПЛІС

Всім привіт!

Одного разу мене запитали замовники, чи немає у мене в проектах цілочисельного ШПФ (Швидке перетворення Фур’є), на що я завжди відповідав, що це вже зроблено іншими у вигляді готових, хоч і кривих, але безкоштовних IP-ядер (Altera / Xilinx) – беріть і користуйтеся. Однак, ці ядра не оптимальні, володіють набором «рис» і потребують подальшої доробки. У зв’язку з чим, пішовши в черговий плановий відпустку, яку не хотілося провести бездарно, я зайнявся реалізацією конфигурируемого ядра цілочисельного ШПФ (Швидке перетворення Фур’є).

КДПВ (процес отдладки помилки переповнення даних)

У статті я хочу розповісти, якими способами і засобами реалізуються математичні операції при обчисленні швидкого перетворення Фур’є в цілочисельному форматі на сучасних кристалах ПЛІС. Основу будь-якого ШПФ представляє вузол, який носить назву «метелик». У метелику реалізуються математичні дії – додавання, множення і віднімання. Саме про реалізацію «метелики» і її закінчених вузлів буде йти розповідь в першу чергу. За основу взяті сучасні сімейства ПЛІС фірми Xilinx – це серія Ultrascale і Ultrascale+, а також зачіпаються старші серії 6- (Virtex) і 7- (Artix, Kintex, Virtex). Більш старші серії в сучасних проектах – не представляють інтересу в 2018 році. Мета статті – розкрити труднощі та особливості реалізації кастомних ядер цифрової обробки сигналів на прикладі ШПФ (Швидке перетворення Фур’є).

Введення

Ні для кого не секрет, що алгоритми взяття ШПФ міцно увійшли в життя інженерів цифрової обробки сигналів, а отже, цей інструмент потрібен постійно. У топових виробників FPGA, таких як Altera / Xilinx вже є гнучкі конфігуровані ядра FFT / IFFT, однак вони мають ряд обмежень і особливостей, у зв’язку з чим мені вже не раз доводилося користуватися власними напрацюваннями. Ось і цього разу мені довелося реалізувати ШПФ в цілочисельному форматі за схемою Radix-2 на ПЛІС. У своїй попередній статті я вже робив ШПФ у форматі з плаваючою точкою, і звідти ви знаєте, що для реалізації ШПФ застосовується алгоритм з дворазовим паралелізмом, тобто ядро може обробляти два комплексних відліки на одній частоті. Це ключова особливість ШПФ, який відсутній у готових ядрах FFT Xilinx.

Приклад: потрібно розробити сайт ШПФ, який здійснює безперервну роботу вхідного потоку комплексних чисел на частоті 800 МГц. Ядро від Xilinx таке не потягне (реальні тактові частоти обробки в сучасних ПЛІС близько 300-400 МГц), або потребуватиме якимось чином децимірувати вхідний потік. Кастомне ядро дозволяє без попереднього втручання тактувати два вхідних відліку на частоті 400 МГц замість одного відліку на 800 МГц. Інший мінус ядра FFT Xilinx пов’язаний з неможливістю приймати вхідний потік в біт-реверсному порядку. У зв’язку з чим витрачається величезний ресурс пам’яті кристала ПЛІС для перестановки даних у нормальний порядок. Для завдань швидкої згортки сигналів, коли два вузла ШПФ стоять один за одним, — це може стати критичним моментом, тобто завдання просто не ляже в обраний кристал ПЛІС. Кастомное ядро ШПФ дозволяє приймати на вході дані в нормальному порядку, а видавати – біт-реверсному, а ядро зворотного ШПФ – навпаки, отримує дані в біт-реверсному порядку, а видає в нормальному. Економиться відразу два буфера на перестановку даних!!!
Оскільки більша частина матеріалу цієї статті могла перетинатися з попередньою, я вирішив викинути все непотрібне, зробивши великий упор на математичні операції в цілочисельному форматі на FPGA для реалізації ШПФ.

Параметри ядра ШПФ

 

  • NFFT – кількість метеликів (довжина ШПФ),
  • DATA_WIDTH – розрядність вхідних даних (4-32),
  • TWDL_WIDTH – розрядність повертаючих множників (8-27).
  • SERIES – задає сімейство ПЛІС, на якій реалізується ШПФ (“NEW” – Ultrascale, “OLD” – 6/7 серія ПЛІС Xilinx).


Як і будь-які інші ланки в ланцюзі, БПФ має вхідні порти управління – тактовий сигнал скидання, а також вхідні та вихідні порти даних. Крім того, в ядрі використовується сигнал USE_FLY, який дозволяє динамічно вимикати метелики ШПФ для процесів налагодження або возжности подивитися вихідний вхідний потік.

У таблиці нижче наведено обсяг займаних ресурсів FPGA в залежності від довжини ШПФ NFFT для DATA_WIDTH = 16 і двох разрядностей TWDL_WIDTH = 16 і 24 біт.

Ядро при NFFT = 64K стабільно працює на частоті обробки FREQ = 375 MHz на кристалі Kintex-7 (410T).

Структура проекту

Для зручності розуміння особливостей тих чи інших компонентів, наведу список файлів проекту і їх короткий опис в ієрархічному порядку:

  • Ядра ШПФ:
    • int_fftNk – вузол ШПФ, схема Radix-2, децимация по частоті (DIF), вхідний потік – нормальний, вихідний – біт-реверсний.
    • int_ifftNk – вузол ОБПФ, схема Radix-2, децімація по часу (DIT), вхідний потік – біт-реверсний, вихідний – нормальний.
  • Метелики:
    • int_dif2_fly – метелик Radix-2, децимация по частоті,
    • int_dit2_fly – метелик Radix-2, децимация за часом,
  • Комплексні помножувачі:
    • int_cmult_dsp48 – загальний конфігурується помножувач, включає в себе:
    • int_cmult18x25_dsp48 – помножувач для невеликих разрядностей даних і повертаючих множників,
    • int_cmult_dbl18_dsp48 – подвоєний помножувач, розрядність повертаючих множників до 18 біт,
    • int_cmult_dbl35_dsp48 – подвоєний помножувач, розрядність повертаючих множників до 25* біт,
    • int_cmult_trpl18_dsp48 – потроєний помножувач, розрядність повертаючих множників до 18 біт,
    • int_cmult_trpl52_dsp48 – потроєний помножувач, розрядність повертаючих множників до 25* біт,
  • Помножувачі:
    • mlt42x18_dsp48e1 – помножувач з разрядностями операндів до 42 і 18 біт на базі DSP48E1,
    • mlt59x18_dsp48e1 – помножувач з разрядностями операндів до 59 і 18 біт на базі DSP48E1,
    • mlt35x25_dsp48e1– помножувач з разрядностями операндів до 35 і 25 біт на базі DSP48E1,
    • mlt52x25_dsp48e1– помножувач з разрядностями операндів до 52 та 25 біт на базі DSP48E1,
    • mlt44x18_dsp48e2 – помножувач з разрядностями операндів до 44 і 18 біт на базі DSP48E2,
    • mlt61x18_dsp48e2 – помножувач з разрядностями операндів до 61 та 18 біт на базі DSP48E2,
    • mlt35x27_dsp48e2 – помножувач з разрядностями операндів до 35 і 27 біт на базі DSP48E2,
    • mlt52x27_dsp48e2– помножувач з разрядностями операндів до 52 і 27 біт на базі DSP48E2.
  • Суматор:
    • int_addsub_dsp48 – універсальний суматор, розрядності операндів до 96 біт.
  • Лінії затримки:
    • int_delay_line – базова лінія затримки, забезпечує перестановку даних між метеликами,
    • int_align_fft – вирівнювання вхідних даних і повертаючих множників на вході метелики ШПФ,
    • int_align_fft – вирівнювання вхідних даних і повертаючих множників на вході метелики ОБПФ,
  • Повертаючи множників:
    • rom_twiddle_int – генератор повертаючих множників, певної довжини ШПФ вважає коефіцієнти на базі DSP осередків ПЛІС,
    • row_twiddle_tay – генератор повертаючих множників за допомогою ряду Тейлора (NFFT > 2K)**.
  • Буфера даних:
    • inbuf_half_path – вхідний буфер, приймає потік в звичайному режимі і видає дві послідовності відліків, зрушених на половину довжини ШПФ***,
    • outbuf_half_path – вихідний буфер, збирає дві послідовності і видає одну безперервну рівну довжині ШПФ,
    • iobuf_flow_int2 – буфер працює в двох режимах: приймає потік в режимі Interleave-2 і видає дві послідовності зрушених на половину довжини ШПФ. Або навпаки, в залежності від опції «BITREV».
    • int_bitrev_ord – простий перетворювач даних з натурального порядку в біт-реверсний.

* — для DSP48E1: 25 біт, для DSP48E2 – 27 біт.
** — з певної стадії ШПФ можна використовувати фіксовану кількість блокової пам’яті для зберігання повертаючих коефіцієнтів, а проміжні коефіцієнти вираховувати з допомогою DSP48 вузлів за формулою Тейлора до першої похідної. У зв’язку з тим, що для ШПФ ресурс пам’яті важливіше, можна сміливо жертвувати обчислювальними блоками на догоду пам’яті.
*** — вхідний буфер і лінії затримки — вносять істотний внесок в обсяг займаних ресурсів пам’яті ПЛІС

«Метелик»
Всі, хто хоча б раз стикалися з алгоритмом швидкого перетворення Фур’є, знають, що в основі цього алгоритму лежить елементарна операція – «метелик». Він перетворює вхідний потік, множачи вхідні дані на повертаючі коефіцієнти (twiddle factor). Для ШПФ є дві класичні схеми перетворення – децимація по частоті (DIF, decimation-in-frequency) і децимація по часу (DIT, decimation-in-time). Для DIT алгоритму характерне розбиття вхідної послідовності на дві послідовності половинної тривалості, а для DIF алгоритму – на дві послідовності парних і непарних відліків тривалістю NFFT. Крім того, ці алгоритми відрізняються математичними діями для операції «метелик».

A, B – вхідні пари комплексних відліків,
X, Y – вихідні пари комплексних відліків,
W – комплексні повертають множники.

Оскільки вхідні дані – комплексні величини, то метелик вимагає один комплексний помножувач (4 операції множення і 2 операції додавання) і два комплексних суматора (4 операції додавання). Це і є весь математичний базис, який необхідно реалізувати на ПЛІС.

Помножувач

Слід зазначити, що всі математичні операції на FPGA часто виконуються в додатковому коді (2’s complement). Помножувач на ПЛІС можна реалізувати двома способами — на логіці, використовуючи тригери і таблиці LUT, або на спеціальних блоках обчислення DSP48, які давно і міцно увійшли до складу всіх сучасних ПЛІС. На логічних блоках множення реалізується з допомогою алгоритму Бута або його модифікацій, займає досить пристойний обсяг логічних ресурсів і далеко не завжди задовольняє тимчасовим обмеженням на високих частотах обробки даних. У зв’язку з цим, помножувачі на ПЛІС в сучасних проектах практично завжди проектуються на базі DSP48 вузлів і лише зрідка – на логіці. Вузол DSP48 – це складна закінчена комірка, яка реалізує математичні та логічні функції. Основні операції: множення, додавання, віднімання, накопичення, лічильник, логічні операції (XOR, NAND, AND, OR, NOR), зведення в квадрат, порівняння чисел, зрушення і т. д. На наступному малюнку представлена осередок DSP48Е2 для сімейства ПЛІС Xilinx Ultrascale+.

Шляхом простої конфігурації портів, операцій обчислення у вузлах і виставлення затримок всередині вузла – можна зробити швидкісну математичну молотарку даних.
Зазначимо, у всіх топових вендорів FPGA в середовищі проектування є стандартні і безкоштовні IP-ядра для обчислення математичних функцій на базі вузла DSP48. Вони дозволяють обчислювати примітивні математичекие функції і виставляти різні затримки на вході і виході вузла. У Xilinx це IP-Core «multiplier» (ver. 12.0, 2018), яке дозволяє конфігурувати помножувач на будь-яку розрядність вхідних даних від 2 до 64 біт. Крім того, можна задати спосіб реалізації помножувача – на логічних ресурсах або вбудованих примітивах DSP48. Оцініть, скільки логіки «їсть» помножувач з розрядністю вхідних даних на портах A та B = 64 біт. Якщо використовувати вузли DSP48, то їх потрібно всього 16.

Основне обмеження, що накладається на клітинки DSP48 — розрядність вхідних даних. У вузлі DSP48E1, який є базовою осередком ПЛІС Xilinx 6 і 7 серії розрядність портів входу для множення: «A» — 25 біт, «B» — 18 – біт, Отже, результат множення являє 43-бітове число. Для сімейства ПЛІС Xilinx Ultrascale і Ultrascale+ вузол зазнав кілька змін, зокрема розрядність першого порту зросла на два біта: «A» — 27 біт, «B» — 18 — біт. Крім того, сам вузол називається DSP48E2.
Щоб не мати прив’язку до конкретного сімейства і кристалу FPGA, забезпечити «чистоту вихідного коду», і врахувати всі можливі розрядності вхідних даних, вирішено спроектувати власний набір умножителей. Це дозволить максимально ефективно реалізувати комплексні помножувачі для «метеликів» ШПФ, а саме — помножувачі і суматор-вычитатель на базі DSP48 блоків. Перший вхід помножувача — це вхідні дані, другий вхід помножувача — повертаючі множники (гармонічний сигнал з пам’яті). Набір умножителей реалізується за допомогою вбудованої бібліотеки UNISIM, з якої необхідно підключити примітиви DSP48E1 і DSP48E2 для подальшого їх використання у проекті. Набір множників представлений в таблиці. Слід зазначити, що:

  • Операція множення чисел призводить до зростання розрядності твори як сума разрядностей операндів.
  • Кожен з множників 25×18 і 27×18 дублюються, по суті – це один компонент для різних сімейств.
  • Чим більше стадія паралельності операції – тим більше затримка на обчислення і більше обсяг позичених ресурсів.
  • При меншої розрядності на вході «B» можна реалізувати помножувачі з більшою розрядністю на іншому вході.
  • Основне обмеження в нарощуванні розрядності вносить порт «B» (реальний порт примітиву DSP48) і внутрішній сдвиговый регістр на 17-біт.

Подальше збільшення розрядності не представляє інтересу в рамках поставленої задачі з причин, описаних нижче:

Розрядність повертаючих множників

Очевидно, що чим більше розрядність гармонічного сигналу – тим точніше представляється число (більше знаків в дробовій частині). Але рзрядность порту B < 25 біт у зв’язку з тим, що для повертаючих множників у вузлах ШПФ цієї розрядності цілком достатньо для забезпечення якісного множення вхідного потоку з елементами гармонічного сигналу в «метеликів» (для будь-яких реально досяжних довжин ШПФ на сучасних ПЛІС). Типове значення розрядності повертаючих коефіцієнтів реалізуються мною завдання 16 біт, 24 – рідше, 32 – ніколи.

Розрядність вхідних відліків

Розрядність даних типових вузлів прийому та реєстрації (АЦП, ЦАП) не велика – від 8 до 16 біт, і рідко – 24 або 32 біта. Причому, в останньому випадку ефективніше скористатися форматом даних з плаваючою точкою за стандартом IEEE-754. З іншого боку, кожна стадія «метелики» ШПФ додає один розряд даних до вихідних відлікам з-за виконання математичних операцій. Наприклад, для довжини NFFT = 1024 використовується log2(NFFT) = 10 метеликів. Отже, вихідна розрядність буде на десять біт більше вхідний, WOUT = WIN + 10. У загальному вигляді формула виглядає так:
WOUT = WIN + log2(NFFT);

Приклад:
Розрядність вхідного потоку WIN = 32 біт,
Розрядність повертаючих множників TWD = 27,
Розрядність порту «А» зі списку реалізованих умножителей в даній статті не перевищує 52 біта. Це означає, що максимальна кількість метеликів ШПФ = 52-32 = 20. Тобто можливо реалізувати ШПФ довжиною до 2^20 = 1М відліків. (Проте, на практиці, прямими засобами таке неможливо у зв’язку з обмеженістю ресурсів навіть у самих потужних кристалів ПЛІС, але це відноситься до іншої теми і в статті розглядатися не буде).

Як видно, це одна з головних причин, по якій я не став реалізовувати помножувачі з більшою розрядністю вхідних портів. Використовувані помножувачі покривають повний діапазон необхідних значень розрядності вхідних даних і повертаючих множників для задачі обчислення цілочисельного ШПФ. У всіх інших випадках можна скористатися обчисленням ШПФ у форматі з плаваючою точкою!

Реалізація «широкого» помножувача

На базі простого прикладу множення двох чисел, що не вписуються в розрядність стандартного сайту DSP48, я покажу, як можна реалізувати широкий помножувач даних. На наступному малюнку представлена його структурна схема. Помножувач реалізує множення двох знакових чисел у додатковому коді, розрядність першого операнда X – 42 біта, другого Y – 18 біт. Вона містить два вузла DSP48E2. Для вирівнювання затримок у верхньому вузлі використовується два регістри. Це зроблено з тієї причини, що у верхньому суматорі потрібно правильно скласти числа з верхнього і нижнього вузлів DSP48. Нижній суматор фактично не використовується. На виході нижнього вузла варто додаткова затримка  для вирівнювання вихідного числа по часу. Сумарна затримка – 4 такти.


Твір складається з двох компонент:

  • Молодша частина: P1 = ‘0’ & X[16:0] * Y[17:0];
  • Старша частина: P2 = X[42:17] * Y[17:0] + (P1 >> 17);

 

Суматор

Аналогічно множнику, суматор може будуватися на логічних ресурсах, використовуючи ланцюжок перенесення, або на DSP48 блоках. Для досягнення максимальної пропускної здатності другий метод краще. Один примітив DSP48 дозволяє реалізувати операцію додавання до 48 біт, два вузла – до 96 біт. Для реалізується завдання таких разрядностей цілком достатньо. Крім того, примітив DSP48 володіє спеціальним режимом «SIMD MODE», який розпаралелює вбудований 48-бітний АЛУ на кілька операцій з різними даними невеликої розрядності. Тобто, в режимі «ONE» використовується повна розрядна сітка в 48 біт і два операнда, а в режимі «TWO» розрядна сітка розбивається на декілька паралельних потоків по 24 біта кожен (4 операнда). Цей режим, використовуючи всього один суматор, допомагає скоротити кількість займаних ресурсів кристала ПЛІС на невеликих разрядностях вхідних відліків (на перших стадіях обчислення).

Зростання розрядності даних

Операція множення двох чисел з разрядностями N і M у двійковому додатковому коді призводить до зростання вихідної розрядності до Р = N + M.
Приклад: для перемноження трехбитовых чисел N = M = 3, максимальне позитивне число дорівнює +3 = (011)2, а максимально негативне число – 4 = (100)2. Старший біт відповідає за знак числа. Отже, максимально можливе число при множенні – це +16 = (010000)2, що утворюється в результаті множення двох максимально від’ємних чисел -4. Розрядність вихідних даних дорівнює сумі вхідних разрядностей Р = N+M = 6 біт.

Операція додавання двох чисел з разрядностями N і M у двійковому додатковому коді призводить до зростання вихідний розрядності на один біт.
Приклад: складемо два позитивних числа, N = M = 3, максимальне позитивне число дорівнює 3 = (011)2, а максимальне від’ємне число – 4 = (100)2. Старший біт відповідає за знак числа. Отже, максимальне додатне число – це 6 = (0110)2, а максимальне від’ємне число – це -8 = (1000)2. Розрядність вихідних даних збільшується на один біт.

Облік особливостей алгоритму

Усічення розрядності зверху:
Для мінімізації ресурсів ПЛІС в алгоритмах ШПФ вирішено при множенні даних у метелику ніколи не використовувати максимально можливе від’ємне число для повертаючих коефіцієнтів. Ця поправка не вносить ніякого негативного впливу на результат. Наприклад, при ипсользовании 16-бітного представлення повертаючих множників мінімальне число -32768 = 0x8000, а наступне за ним -32767 = 0x8001. Похибка при заміні максимально від’ємного числа найближчим сусіднє значення складе ~0.003% і повністю прийнятним для реалізації завдання.
Прибравши в добутку двох чисел мінімальне число, на кожній ітерації можна скоротити один невикористаний старший розряд. Приклад: дані – 4 = (100)2, коефіцієнт +3 = (011)2. Результат множення = -12 = (110100)2. П’ятий біт можна відкинути, оскільки він дублює сусідній, четвертий — знаковий біт.

Усічення розрядності знизу:
Очевидно, що примножуючи в «метелику» вхідний сигнал на гармонійний вплив, не потрібно тягнути вихідну розрядність у наступні метелики, а потрібно округлення або скорочено. Повертаючи множників представлені в зручному M-бітному форматі, однак насправді це звичайний синус і косинус, нормовані до одиниці. Тобто число 0x8000 = -1, а число 0x7FFF = +1. Отже, результат множення обов’язково скорочується до вихідної розрядності даних (тобто M біт від повертаючих множників скорочуються знизу). У всіх реалізаціях ШПФ, що мені довелося бачити, повертаючи множників нормовані до 1 тими чи іншими способами. Отже, з результату множення необхідно взяти біти в такій сітці [N+M-1-1: M-1]. Старший біт не використовується (віднімемо додаткову 1), молодші – скорочуються.

Додавання/віднімання даних в операціях «метелики» ніяк не піддається мінімізації і тільки ця операція вносить внесок у зростання розрядності вихідних даних на один біт на кожній стадії обчислення.

Зазначимо, що на першій стадії алгоритму FFT DIT або на останній стадії алгоритму FFT DIF дані необхідно помножити на повертаючий множник з нульовим індексом W0 = {Re, Im} = {1, 0}. У зв’язку з тим, що множення на одиницю і нуль – примітивні операції, їх можна не виконувати. В такому випадку операція комплексного множення взагалі не потрібно: реальна і уявна компоненти проходять «поворот» без змін. На другій стадії використовується два коефіцієнта: W0 = {Re, Im} = {1, 0} і W1 = {Re, Im} = {0, -1}. Аналогічно, операції можна звести до елементарних перетворень і використовувати мультплексор для вибору вихідного відліку. Це дозволяє істотно заощадити DSP48 блоки на перших двох метеликів.

Комплексний помножувач будується аналогічно – на базі множників і суматора-від’ємника, однак для деяких варіантів розрядності вхідних даних додаткові ресурси не потрібні, про що буде розказано нижче.

Вхідний буфер і лінії затримки і крос-комутатори аналогічні тим, що описані в попередній статті. Повертаючи множників стали цілочисельними з конфігурованою розрядністю. По іншому — глобальних змін у проектуванні ядра ШПФ — нема

Особливості ядра ШПФ INT_FFTK

 

  • Повністю конвеєрна схема обробки даних.
  • Довжина перетворення NFFT = 8-512K точок.
  • Гнучка настройка довжини перетворення NFFT.
  • Цілочисельний формат вхідних даних, розрядність конфігурується.
  • Цілочисельний формат повертаючих множників, розрядність конфігурується.
  • Компактне зберігання повертаючих коефіцієнтів через розкладання в ряд Тейлора до першої похідної.
  • Зростання розрядності на кожній стадії обчислення і на виході ядра!
  • Зберігання чверті періоду коефіцієнтів для економії ресурсів кристала.
  • ШПФ: дані на вході – у прямому порядку, на виході в двійково-інверсному.
  • ОБПФ: дані на вході в двійково-інверсному порядку, на виході – у прямому.
  • Конвеєрна схема обробки з дворазовим паралелізмом. Radix-2.
  • Мінімальна нерозривний довжина пачки даних дорівнює NFFT відліків*.
  • Висока швидкість обчислень і невеликий обсяг позичених ресурсів.
  • Реалізація на останніх кристалах ПЛІС (Virtex-6, 7-Series, Ultrascale).
  • Частота роботи ядра ~375MHz на Kintex-7
  • Мова опису – VHDL.
  • Відсутність необхідності операції перетворення bitreverse для зв’язки ШПФ+ОБПФ.
  • Open Source проект без включення сторонніх IP-Cores.

 

Вихідний код

Вихідний код ядра ШПФ INTFFTK на VHDL (включаючи базові операції і набір умножителей) і m-скрипти для Matlab / Octave доступні в моєму профілі на гітхабі.

Висновок

В ході розробки спроектовано нове ядро ШПФ, що забезпечує високу продуктивність у порівнянні з аналогами. Зв’язка ядер ШПФ і ОБПФ не вимагає перекладу в натуральний порядок, а максимальна довжина перетворення обмежена лише ресурсами ПЛІС. Дворазовий паралелізм дозволяє обробляти вхідні потоки подвоєної частоти, чого не може робити IP-CORE Xilinx. Розрядність на виході целочисленого ШПФ лінійно зростає у залежності від кількості стадій перетворення.
Спасибі за увагу!

Related Articles

Залишити відповідь

Ваша e-mail адреса не оприлюднюватиметься. Обов’язкові поля позначені *

Close