Всем привет! В этой статье я хочу рассказать о реализации сверхдлинного алгоритма быстрого преобразования Фурье на FPGA. На написание этой статьи меня побудило желание поделиться личным практическим опытом, который мне не хотелось бы потерять, оставив информацию только в голове.
А поскольку я больше не занимаюсь задачами цифровой обработки сигналов на ПЛИС, мне просто приходится передавать имеющиеся у меня знания.
В данной статье показана невозможность реализации «классической» очень длинной схемы БПФ даже на самых современных микросхемах FPGA и предложен алгоритм, позволяющий это сделать.
Также поэтапно рассматривается основная идея алгоритма: от математической составляющей до создания полноценного решения на базе ПЛИС с использованием внешней памяти DDR. В статье будут затронуты тонкости проектирования многоканальных систем обработки для данного класса задач и, в частности, описан мой практический опыт. 
Постановка задачи
Представьте, что вам нужно спроектировать широкополосный анализатор спектра с высоким частотным разрешением.Известно, что чем больше длина БПФ, тем выше разрешение.
То есть нужно спроектировать систему, которая принимает сигналы от АЦП, выполняет некоторую обработку внутри ПЛИС и выводит данные через высокоскоростной интерфейс (например, PCIe, USB и т.д.) на следующие этапы обработки данных.
.
Формализуем требования задачи:
- Алгоритм: быстрое преобразование Фурье
- Длина БПФ от 256 тыс.
до 256 млн точек.
- Непрерывная обработка данных в режиме реального времени.
- Количество независимых каналов на входе от 1 до N.
- Ширина сигнала АЦП от 16 бит.
- Комплексный сигнал на входе.
- Частота данных, принимаемых с одного канала АЦП, от 200 МГц и выше.
- Наличие переключаемой оконной функции перед звеном БПФ.
Классический подход не работает!
Прежде чем перейти к описанию алгоритма сверхдлинного преобразования Фурье, следует сказать о принципиальной невозможности реализации классических схем БПФ определенной длины на ПЛИС.Именно из-за этого и возникает решение, которое будет описано в этой статье.
Представим, что перед вами стоит задача разработать узел БПФ длиной от 256К до 64М отсчетов.
В свободном доступе нет ядер с длиной БПФ, превышающей 64 тыс.
точек, от ведущих производителей Xilinx и (например, Altera) Intel. Причина – огромное потребление блочной памяти (BRAM) кристалла ПЛИС при увеличении длины БПФ.
Не помогает и появление блоков URAM в микросхемах от Xilinx. Даже если вы напишете свое сверхоптимальное решение с использованием классического алгоритма БПФ через конвейерное соединение узлов-бабочек и узлов кросс-связи, реализовать его на ПЛИС у вас не получится.
Почему? Когда длина БПФ увеличивается в N раз, количество занятых ресурсов BRAM изменяется пропорционально, по крайней мере, на такую же величину.
Для БПФ с плавающей запятой количество узлов BRAM увеличивается пропорционально, но для БПФ с фиксированной запятой оценка потребления ресурсов еще более пессимистична.
Это связано с тем, что на каждом этапе расчета БПФ (после каждой бабочки) разрядность промежуточных данных увеличивается на 1 бит.
Итоговая разрядность на выходе узла БПФ рассчитывается по формуле 
Например, БПФ длиной 1 млн выборок добавляет 20 бит к входному сигналу.
Если входной сигнал имеет ширину 16 бит, то выходной сигнал будет соответственно 36 бит. 
Давайте посмотрим на таблицу занятых ресурсов для выборок FFT 64K на примере Xilinx. По приблизительным оценкам, если вам нужно создать блок БПФ с плавающей запятой из 1 млн точек, вам понадобится 400 * 16 = 6400 элементы блочной памяти типа RAMB36K или ~ 220 Мбит ! Такими ресурсами обладает только топовый Virtex UltraScale+ (VU29P), и то за счет ячеек URAM. 
В следующих таблицах показано, сколько внутренней памяти имеют современные ПЛИС от Xilinx, на примере этой серии.
Кинтекс .
Кинтекс Ультрамасштаб 
Кинтекс Ультрамасштаб+ 
Как видите, ресурсы ПЛИС для блочной памяти весьма ограничены, поэтому реализовать БПФ даже для длины 1М точек в классической форме практически невозможно.
невозможный .
Очень длинное БПФ
К счастью, вы можете использовать небольшой математический трюк, чтобы получить одномерное БПФ из двумерного.Общая идея алгоритма заключается в том, что вектор сигнала длиной N разбивается на отсчеты N1 и N2 (где N1 и N2 кратны степеням двойки).
Этот вектор преобразуется в матрицу размерами N1 x N2, над которой производятся все расчеты.
Короткие БПФ длиной N1 и N2 применяются к строкам и столбцам.
Формула расчета БПФ: 
Разделите последовательность N на произведение N1 и N2, затем: 
Тогда формула примет вид: 
После преобразования поворотных коэффициентов: 
Результат: 
Вы можете видеть, что одномерное БПФ длины N превращается в два БПФ длиной N1 и N2, а результат первого БПФ умножается на вращающиеся коэффициенты.
Блок-схема алгоритма сверхдлинного БПФ: 
Этапы алгоритма:
- Транспонирование одномерной последовательности в матрицу K1 x K2
- Вычислить БПФ по K1 строкам K2 раз
- Транспонирование полученной матрицы
- Умножение результатов на поворотные коэффициенты
- Вычислить БПФ по K2 столбцам K1 раз
- Транспонирование матрицы в исходный одномерный вектор
Пример: Ядро БПФ длиной 1М выборок.
Вы можете разделить вычисление на два выборки БПФ по 1 КБ: 1024 x 1024 = 1048576. На следующем рисунке показано, что для узла БПФ с 1024 точками требуется только 7 ячеек RAMB36K. 
Видно, что небольшие ядра БПФ практически не занимают ресурсы ПЛИС, в частности, почти не используют блочную память.
Давайте пройдемся по всем элементам алгоритма и рассмотрим основные особенности каждого.
Пошаговая реализация
Блоки БПФ
Узлы БПФ реализованы по классической конвейерной схеме – последовательное соединение log2(N) бабочек Radix-2/Radix-4 и узлов кросс-соединения.Вы можете использовать готовые бесплатные ядра от производителей FPGA или написать собственное оптимизированное ядро.
Узлы БПФ могут выполнять вычисления в формате с фиксированной, плавающей запятой или в каком-либо другом формате, в зависимости от задачи.
Узел БПФ:
- Конвейерный, потоковый ввод-вывод
- Фиксированная немасштабированная/плавающая точка
- Радикс-2 / Радикс-4
- Обратный бит/естественный порядок
- Комплексный множитель: структура 4-DSP
Возникает вопрос: устраивает ли вас точность вычислений выбранного формата данных или нет? В моей практике для сверхдлинных БПФ от длин 1М ошибки округления операций с плавающей запятой (IEEE 754) портили результирующий спектр из-за маленькой мантиссы (25 бит).
Поэтому необходимо было использовать БПФ в формате с фиксированной запятой на всех этапах обработки.
В частном случае можно написать свою собственную математику с плавающей запятой на блоках FPGA DSP48 с расширенной мантиссой (результат того стоит).
Ядро БПФ может быть реализовано для вывода результата в обратном порядке битов или в естественном (нормальном) порядке.
В первом случае вы экономите еще больше ресурсов блочной памяти, но немного усложняете алгоритмы перестановки данных при транспонировании матриц.
Если вы только начали реализовывать сверхдлинные БПФ, начните со второго варианта, потому что его легче отлаживать.
Если вы реализуете многоканальную систему, и ядро БПФ вы написали самостоятельно, не используя готовое решение от производителей вашей ПЛИС, то вы можете дополнительно сэкономить память BRAM на хранении коэффициентов поворота для бабочек БПФ.
Для N независимых параллельных каналов можно использовать только 1 модуль хранения с вращающимся умножителем.
То есть, чем больше системный канал, тем больше экономия ресурсов.
Вращающиеся множители
Согласно схеме алгоритма, перед подачей матрицы во второе звено БПФ данные необходимо умножить на коэффициенты поворота.Есть два способа сделать это: использовать DDS или воспользоваться алгоритмом CORDIC. Первый способ предполагает хранение большого объема данных, что требует значительного объема блочной памяти ПЛИС, за что мы бьемся с самого начала этой статьи.
Теоретически можно использовать приближение Тейлора и уменьшить хранимый массив в BRAM. Но в моей практике такое решение искажает результирующий спектр из-за ступенчатой формы сигнала поворотных умножителей.
Второй метод на основе CORDIC вообще не требует блочной BRAM, поскольку использует итеративную схему применения операций сдвига и сложения/вычитания.
К недостаткам алгоритма CORDIC можно отнести длительное время расчета следующего значения (требуется около 20-30 тактов, количество зависит от разрядности).
Этот недостаток приводит к организации дополнительной линии задержки входящих данных, что отнимает определенный логический ресурс ПЛИС.
Например, для многоканальной схемы разрядностью 512 (2 комплексных выборки по 32 бита, 8 каналов) дополнительно потребуется 512*30 = 15 тысяч триггеров .
Для этого в ПЛИС Xilinx имеются ячейки SLICEM, организующие линию задержки на сдвиговых регистрах.
Либо можно потратить на линию задержки несколько блоков BRAM. К ядру CORDIC выдвигаются следующие требования:
- Схема параллельной обработки
- Разрядность входных данных должна быть не менее log2(N) конечной длины БПФ, чтобы обеспечить отсутствие искажений в спектре сигнала.
Для многоканальной схемы также можно использовать один поворотный узел умножителя на весь поток данных и сэкономить ресурсы ПЛИС.
Транспонировать узлы
Контроллеры памяти используются для хранения векторов промежуточных данных, а также для транспонирования матрицы на всех этапах алгоритма.Это может быть любая доступная вам память: QDR SRAM, DDR3 или DDR4 SDRAM. Последние два я использовал в своих проектах.
Но общие принципы работы те же: контроллер памяти транспонирует выборку — получает пакет данных» построчно ", но создает пакет данных в формате " по столбцам ".
Как видно из схемы алгоритма, для данной задачи необходимы три узла внешней памяти, к которым предъявляются два основных требования.
Первый: Память должна хранить минимум выборку длиной в 2 длины БПФ.
Это необходимо для того, чтобы иметь возможность дочитать второй пакет в обратном виде (по столбцам матрицы) при записи одного пакета данных в прямом виде (по строкам матрицы).
Но лучшее решение для хранения данных — не через мультиплексор, а при реализации внешней памяти.
по схеме ФИФО .
В этом случае внешняя память может хранить множество пакетов данных длиной N и эффективно использовать свой ресурс.
Также на практике такая схема позволяет бороться с небольшими замираниями интерфейса на выходе узла БПФ.
В частности, при кратковременном замирании передачи по шине PCIe вероятность переполнения памяти в режиме FIFO существенно ниже, чем в схеме переключения между одним пакетом и другим.
В проектах, которые я реализовал, память DDR почти всегда переполнялась при зависании шины PCIe в режиме мультиплексора, но никогда в режиме FIFO. Рассчитаем объём хранения данных во внешней памяти.
Пусть ширина входных данных составляет 32 бита (одиночное число с плавающей запятой), сигнал комплексный (I/Q), длина БПФ — 1М отсчетов, схема реализации — « настольный теннис " как минимально необходимое требование.
Тогда для хранения потребуется 2*32/8 = 8 МБ Память.
Для хранения данных в режиме FIFO глубиной 32 вам понадобится 256 МБ памяти .
Второй: Память должна иметь возможность записывать и читать данные в реальном времени.
По алгоритму схема передачи данных на входе и выходе контроллера различна.
Поэтому необходимо правильно организовать процесс передачи данных на максимальной скорости, чтобы не было обрывов и повреждения данных.
Для этого еще на этапе проектирования или покупки готовой платы, на которой будет работать сверхдлинное БПФ, нужно рассчитать максимальную пропускную способность всей системы.
То есть найти максимальный поток ввода-вывода во внешнюю память по даташиту и определить скорость подачи данных с многоканального АЦП.
Пропускная способность не зависит от длины БПФ.
Кроме того, чип FPGA необходим для взаимодействия с тремя контроллерами памяти.
Например, для каждого контроллера памяти SO-DIMM DDR4 SDRAM x64 требуется три банка FPGA (что эквивалентно примерно 150 физическим контактам ввода-вывода кристалла).
Всего потребуется не менее 450 портов ввода-вывода или 9 банков HP (High-Performance) FPGA, не считая банков мультигигабитных трансиверов и банка конфигурации Bank0.
Пример настроек IP-ядра DDR4 SDRAM: 
Алгоритм взаимодействия с контроллером памяти:
- входной поток поступает в небольшой FIFO на частоте обработки и преобразуется в поток на частоте контроллера памяти (например, 333 МГц)
- Данные записываются в память в прямой форме, при этом счетчик адреса памяти DDR увеличивается линейно.
- после записи пакета размером N=K1*K2 чтение происходит в обратной форме.
- Счетчик адреса памяти DDR для чтения изменяется по определенному алгоритму с заданным шагом, чтобы обеспечить транспонирование матрицы К1 х К2.
- одновременно с чтением продолжается запись в другую, свободную область внешней памяти.
- поток считываемых данных преобразуется из частоты контроллера в частоту обработки (с использованием небольшого FIFO).
- если по каким-то причинам происходит остановка (замирание) на стороне чтения, запись не прекращается и продолжается до тех пор, пока контроллер не установит флаг ПОЛНЫЙ ФИФО .
Следует отметить, что для достижения максимальной производительности при чтении необходимо пройти все группы FSM контроллера памяти в ПЛИС (см.
даташит на Xilinx DDR MIG), то есть нужно пройти все банки и группы физических банки памяти.
Это накладывает дополнительные ограничения и приводит к необходимости наличия буфера данных после контроллера памяти в ПЛИС.
Его цель — организовать и упаковать обработанные транзакции из каждого банка памяти для дальнейшей передачи на следующие узлы схемы.
Модули URAM (блоки глубиной 4К и шириной 72 бита), появившиеся в современных ПЛИС от Xilinx, идеально подходят для реализации буфера данных.
Также отмечу, что последний узел транспонирования можно использовать как аккумулятор, реализуя схему накопление спектральных составляющих усреднить полученный спектр.
Функция окна
Окно входного сигнала не является обязательным, но если вам нужно добавить оконную функцию к длинному БПФ, вы можете столкнуться с трудностями.Поскольку мы используем сверхдлинные БПФ, нам совершенно негде хранить массив для оконных функций.
БПФ длины N потребует много ресурсов блочной памяти.
Например, оконная функция в виде 32-битного сигнала длиной 16М отсчетов потребует с учетом ее симметрии: 32*4/2=256 Мбит .
Это очень много даже для топовых кристаллов FPGA. Что делать, если вам необходимо иметь возможность непрерывного переключения функций (потребуются как минимум два независимых буфера данных)? Вы можете решить эту проблему очень просто, используя Окна Блэкмана-Харриса необходимый порядок и стандартная формула расчета оконных коэффициентов.
Используя известный CORDIC для генерации гармонических сигналов нужной частоты, можно реализовать оконную функцию Блэкмана-Харриса любого порядка без использования блочной памяти ПЛИС! 
Чем выше порядок оконной функции, тем ниже уровень боковых лепестков спектра.
На практике мне приходилось использовать окна до 5 порядка.
Не будем на этом останавливаться, подробнее о фильтрации окон на ПЛИС.
Я уже говорил об этом в своей предыдущей статье.
Контрольно-пропускные пункты
Поток сигнала через узлы алгоритма сверхдлинного БПФ показан ниже.В качестве входного воздействия выбран сигнал в виде пика при одном значении и ЛЧМ-сигнал.
«Дельта-функция»: 
Чирп-сигнал: 
Скрипт, генерирующий данные на каждом этапе сверхдлинного БПФ.
Написанный на Python, он требует наличия библиотек numpy, scipy и matplotlib. Код Python для изображений
from datetime import datetime from functools import lru_cache import matplotlib.pyplot as plt import numpy as np from scipy.fftpack import fft N1 = 8 # Colomn (FFT1) N2 = 16 # Rows (FFT2) class SignalGenerator: def __init__( self, nfft: int, freq: float, alpha: float = 0.01, ): """Generate some useful kind of signals: harmonic or linear freq. modulated. Parameters ---------- nfft : int Total lenght of FFT (NFFT = N1 * N2).freq : float Signal frequency. alpha : float Add Gaussian noise if alpha not equal to zero. Should be positive. """ self.nfft = nfft self.freq = freq self.alpha = alpha def awgn(self): np.random.seed(42) return self.alpha * np.random.rand(self.nfft) def input_harmonic(self): """Generate input singal""" return ( np.array([1 + 1j if i in (self.freq, self.nfft - self.freq) else 0 for i in range(self.nfft)]) + self.awgn() ) def input_linfreq(self): tt = np.linspace(0, 1, self.nfft, endpoint=False) sig_re = np.cos(self.freq * tt ** 2 * np.pi) * np.sin(tt * np.pi) + self.awgn() sig_im = np.sin(self.freq * tt ** 2 * np.pi) * np.sin(tt * np.pi) + self.awgn() return sig_re + 1j * sig_im class UltraLongFFT: """Calculate ultra-long FFT via 2D FFT scheme with 3 shufflers: step by step. Parameters ---------- n1 : int Rows (number of points for 1st FFTs) n2 : int Columns (number of points for 2ns FFTs) where NFFT = N1 * N2 """ _plt_titles = ( "1. Input Signal", "2. Shuffle [1]", "3. FFT1, N1", "4. Shuffle [2]", "5. Twiddles", "6. Complex Multiplier", "7. FFT2, N2", "8. Shuffle [3].
Output", ) def __init__(self, n1: int = 32, n2: int = 32): self.n1 = n1 self.n2 = n2 self.__nfft = n1 * n2 @property @lru_cache(maxsize=4) def twiddles(self): """Twiddle factors for rotating internal sequence.""" twd = np.array( [np.exp(-1j * 2 * np.pi * (k1 * k2) / self.__nfft) for k1 in range(self.n1) for k2 in range(self.n2)] ) return np.reshape(twd, (self.n1, self.n2)) def fft_calculate(self, signal: np.ndarray) -> np.ndarray: """Calculate signals for each stage of Ultra-long FFT Algorithm Parameters ---------- signal : np.ndarray Input signal. Can be complex. Returns ------- np.ndarray List of arrays for each stage of Ultra-long FFT. """ # ---------------- ULFFT Algorithm ---------------- # 1 Step: Shuffle input sequence sh0_data = np.reshape( a=np.array([signal[k2 * self.n1 + k1] for k1 in range(self.n1) for k2 in range(self.n2)]), newshape=(self.n1, self.n2), ) # 2 Step: Calculate FFT N1 and shuffle res_fft0 = np.array([fft(sh0_data[k1, .
]) for k1 in range(self.n1)]) # 3 Step: Complex multiplier cmp_data = res_fft0 * self.twiddles # 4 Step: Calculate FFT N2 and shuffle res_fft1 = np.array([fft(cmp_data[.
, k2]) for k2 in range(self.n2)]) # Internal Sequences: return np.vstack( [ signal, sh0_data.reshape(-1, self.__nfft), res_fft0.reshape(-1, self.__nfft), res_fft0.T.reshape(-1, self.__nfft), self.twiddles.reshape(-1, self.__nfft), cmp_data.T.reshape(-1, self.__nfft), res_fft1.reshape(-1, self.__nfft), res_fft1.T.reshape(-1, self.__nfft), ] ) def plot_result(self, data: np.ndarray, save_fig: bool = False): """Plot signals for each stage of Ultra-long FFT Algorithm""" _ = plt.figure("Ultra-long FFT", figsize=(16, 12)) for i, arr in enumerate(data): plt.subplot(4, 2, i + 1) plt.plot(arr.real, linewidth=1.15, color="C2") plt.plot(arr.imag, linewidth=1.15, color="C4") plt.title(self._plt_titles[i], fontsize=14) plt.grid(True) plt.xlim([0, self.__nfft - 1]) plt.tight_layout() if save_fig: plt.savefig(f"figure_{datetime.now()}"[:-7]) plt.show() if __name__ == "__main__": _input = SignalGenerator(N1 * N2, freq=16, alpha=0.001) # _array = _input.input_harmonic() _array = _input.input_linfreq() _ulfft = UltraLongFFT(N1, N2) _datas = _ulfft.fft_calculate(_array) _ulfft.plot_result(_datas, save_fig=False)
Чарты в Вивадо
Транспонированный ЛЧМ-сигнал до и после одного из этапов БПФ:
Прохождение ЛЧМ-сигнала через узел БПФ:
Тематическое исследование
Предположим, вы столкнулись с задачей со следующими входными данными:- Входной сигнал сложный, 16-битный.
- Количество каналов – 1. Частота дискретизации 400 МГц.
- Промежуточные данные в формате IEEE-754.
- Длина БПФ: 16 миллионов точек.
Нет функции окна.
Промежуточный информационный поток: ~3 ГБ/с.
Объем памяти для хранения одного пакета БПФ: (2*32/8)*16М = 128 МБ .
Пропускная способность памяти на чтение и запись – не менее ~6 ГБ/с .
Этому требованию соответствует DDR4-2400 SDRAM x32 по формуле: Частота * Двойная скорость * (Ширина / Байт) = 1,2e9 * 2 * (32/8) / 2^30 = 8,94 ГБ/с .
Используя IP Catalog, мы создадим ядра DDR, FFT, CORDIC. Полное БПФ из 16 миллионов точек разлагается в матрицу БПФ размером 4 х 4 КБ.
Предположим, что для переупаковки данных до и после контроллера DDR требуются FIFO, занимающие 4 ячейки блочной памяти RAMB36K. Грубая оценка потребления ресурсов ПЛИС для реализации алгоритма.
Нам нужны 3 контроллера памяти, 2 узла БПФ по 4К, один CORDIC, 6 блоков FIFO (до и после контроллеров памяти), 3 буфера банка памяти.
Для простоты мы не будем учитывать линии задержки для согласования потоков, комплексные множители и другую логику.
Итоговая таблица: 
Как видите, проект занимает не очень много ресурсов и размещается на относительно дешевых ПЛИС.
Однако не стоит забывать, что для 3-х контроллеров памяти требуется определенное количество портов ввода-вывода ПЛИС, поэтому не все микросхемы подходят.
Заключение
В данной статье показан метод реализации сверхдлинных узлов БПФ на ПЛИС для задач высокоточного анализа спектра сигналов в реальном времени.Проектирование такого алгоритма требует определенной настройки «аппаратной» части — к ПЛИС необходимо подключить три независимых контроллера памяти.
Однако это позволяет создавать очень длинные схемы БПФ с количеством точек от 256К до 256М.
Поскольку алгоритм имеет множество нюансов в аппаратной реализации, необходимо заранее рассчитать параметры всех узлов схемы и убедиться в работоспособности ядра на выбранной вами конфигурации.
Платы идеально подходят для реализации сверхдлинного алгоритма БПФ.
Альвео У200/У250 от Xilinx (4 контроллера DDR4 на борту) или Альвео U280 (два модуля DIMM DDR4 и HBM2).
Библиография
- Авторский курс лекций по цифровой обработке сигналов
- Особенности фильтрации окон на FPGA
- Айфичер?.
, Джервис Б.
Цифровая обработка сигналов.
Практический подход
- Рабинер и Голд. Теория и практика цифровой обработки сигналов
- Xilinx PG109: быстрое преобразование Фурье
- Xilinx PG105: КОРДИК
- Xilinx PG150: сверхмасштабная сварка MIG
- Сверхдлинное IP-ядро БПФ для ПЛИС Xilinx
- Сверхдлинная архитектура БПФ
-
Просмотр Видео Mkv На Mac Os X Lion
19 Oct, 24 -
Мой Питомец - Linguaplayer
19 Oct, 24 -
Робот-Гуманоид За 300 Долларов
19 Oct, 24 -
Полноценный Интернет-Магазин На Вашем Сайте
19 Oct, 24
