Здесь приведен полный перевод на русский язык части 3 оригинальной статьи "A Software-Defined Radio for the Masses, Part 3" By Gerald Youngblood, AC5OG. Описывается полная программно-определяемая (SDR) радиостанция на базе ПК, использующая звуковую карту и инновационную схему детектора. В объяснении математика сведена к минимуму. Перевод первой части см. в [18]. Незнакомые термины и сокращения см. в Словарике части 1 [18].

В части 1 этой серии статей [18] было предоставлено общее описание цифровой обработки сигналов (digital signal processing, DSP), как она используется в программно-опредеяемых радиоустройствах (software-defined radios, SDRs), включая обзор полнофункционального трансивера, который основан на PC. Все функции DSP и функции управления реализованы в программе на языке Visual Basic. Часть 2 описывает исходный код Visual Basic, в котором реализован полнодуплексный квадратурный интерфейс (full-duplex quadrature interface) со звуковой картой на PC.

Как было описано ранее in-phase (I) сигнал и quadrature (Q) сигнал дают возможность виртуально модулировать или демодулировать сигналы любого типа. Детектор Tayloe, описанный в части 1 [18], представляет собой простой метод преобразования модулированного RF-сигнала квадратурную форму базового диапазона частот звуковой карты, так что сигналы I и Q можно подать на левый и правый каналы stereo звуковой карты PC для дальнейшей программной обработки. Full-duplex интерфейс DirectX8, описанный в части 2 [19], осуществляет ввод и вывод оцифрованных квадратурных сигналов через звуковую карту. Интерфейс звуковой карты обеспечивает входной буферный массив inBuffer() и выходной буферный массив outBuffer(), через который код DSP принимает захваченный сигнал и затем выводит обработанные данные сигнала.

Эта статья расширяет программный интерфейс звуковой карты для демонстрации реализации приемника SDR. Для этого в программном обеспечении реализованы следующие функции:

• Разделение буферов стереозвука на каналы I и Q.
• Преобразование сигнала из домена времени (time domain) в домен частот с помощью быстрого преобразования Фурье (FFT).
• Декартово-полярное (cartesian-to-polar) преобразование векторов сигнала.
• Преобразование частоты из ПЧ основной полосы со смещением 11.25 кГц в 0 Гц.
• Выбор боковой полосы сигнала.
• Генерация коэффициента полосового фильтра.
• Быстрая сверточная фильтрация FFT (fast-convolution filtering).
• Преобразование сигнала обратно в time domain с помощью обратного быстрого преобразования Фурье (IFFT).
• Цифровая автоматическая подстройка усиления (automatic gain control, AGC) с переменным временем реакции.
• Передача обработанного сигнала в выходной буфер для операций радиопередачи или радиоприема.

Демонстрационный исходный код можно загрузить с сайта ARRLWeb [3]. Программа требует динамически линкуемой библиотеки (dynamic link library, DLL) из Intel Signal Processing Library [4], размещенной в рабочей директории проекта и скомпилированной программы. Эти файлы включены в демонстрационную программу.

[Архитектура программного обеспечения]

На рис. 1 приведена блок-схема архитектуры программного обеспечения DSP. Эта архитектура одинаково хорошо работает как для радиоприема, так и для радиопередачи, причем для переключения между этими режимами требуется изменить всего лишь несколько строк кода. Хотя на блок-схеме показаны модули коррекции амплитуды и фазы (Amplitude and Phase Correction) и адаптивного цифрового фильтра для уменьшения шума и режекции нежелательных сигналов (LMS Noise & Notch Filter), обсуждение этих фич выходит за рамки статьи.

Рис. 1. Блок-схема программной архитектуры DSP.

Коррекция амплитуды и фазы позволяет корректировать несовершенства в фазовом и амплитудном дисбалансе, созданном в аналоговой схеме, в частотной области. LMS noise and notch filters [5] являются адаптивной формой фильтрации с конечной импульсной характеристикой (FIR), которая обеспечивает снижение шума во временной области. Существуют другие методы подавления шума, которые могут быть выполнены в частотной области, такие как спектральное вычитание [6], корреляция [7 [и усреднение FFT [8].

[Парсинг входных буферов для получения векторов сигналов I и Q]

Рис. 2 предоставляет декларации переменных и констант демонстрационного кода.

Public Const Fs As Long = 44100                    'Частота дискретизации, выборок в секунду
Public Const NFFT As Long = 4096                   'Количество бинов FFT
Public Const BLKSIZE As Long = 2048                'Количество выборок в блоке захвата/проигрывания (capture/play)
Public Const CAPTURESIZE As Long = 4096            'Количество выборок в Capture Buffer
Public Const FILTERTAPS As Long = 2048             'Количество коэффициентов (taps) в полосовом фильтре
Private BinSize As Single                          'Размер бинов FFT в Гц, т. е. разрешающая способность FFT

Private order As Long                              'Calculate Order power of 2 from NFFT
Private filterM(NFFT) As Double                    'Polar Magnitude of filter freq resp
Private filterP(NFFT) As Double                    'Polar Phase of filter freq resp
Private RealIn(NFFT) As Double                     'Буферы FFT
Private RealOut(NFFT) As Double
Private ImagIn(NFFT) As Double
Private ImagOut(NFFT) As Double

Private IOverlap(NFFT - FILTERTAPS - 1) As Double  'Overlap prev FFT/IFFT
Private QOverlap(NFFT - FILTERTAPS - 1) As Double  'Overlap prev FFT/IFFT

Private RealOut_1(NFFT) As Double                  'Fast Convolution Filter buffers
Private RealOut_2(NFFT) As Double
Private ImagOut_1(NFFT) As Double
Private ImagOut_2(NFFT) As Double

Public FHigh As Long                               'Верхняя частота среза в Гц
Public FLow As Long                                'Нижняя частота среза в Гц
Public Fl As Double                                'Нижняя частота среза, выраженная как доля (коэффициент)
                                                   ' от частоты дискретизации Fs.
Public Fh As Double                                'Верхняя частота среза, выраженная как доля (коэффициент)
                                                   ' от частоты дискретизации Fs.
Public SSB As Boolean                              'True для режимов однополосной модуляции
Public USB As Boolean                              'Переменная для выбора полосы SSB
Public TX As Boolean                               'True для выбора режима передачи
Public IFShift As Boolean                          'True для ПЧ 11.025 кГц

Public AGC As Boolean                              'True, если АРУ разрешена
Public AGCHang As Long                             'AGC AGCHang time factor
Public AGCMode As Long                             'Сохраняет выбор режима AGC
Public RXHang As Long                              'Сохраняет установку времени RX Hang
Public AGCLoop As Long                             'Счетчик буфера времени AGC AGCHang
Private Vpk As Double                              'Пик отфильтрованного выходного сигнала
Private G(24) As Double                            'Буфер времени усиления AGCHang
Private Gain As Double                             'Установка состояния усиления для AGC
Private PrevGain As Double                         'Усиление AGC во время предыдущего входного блока
Private GainStep As Double                         'Шаги времени атаки AGC
Private GainDB As Double                           'Усиление AGC в dB
Private TempOut(BLKSIZE) As Double                 'Временный буфер для вычисления усиления
Public MaxGain As Long                             'Максимальный коэффициент усиления AGC

Private FFTBins As Long                            'Количество бинов FFT для отображения
Private M(NFFT) As Double                          'Полярная магнитуда двойной точности
Private P(NFFT) As Double                          'Угол фазы двойной точности
Private S As Long                                  'Счетчик циклов для выборок

Private order As Long   'Calculate Order power of 2 from NFFT

Private filterM(NFFT) As Double   'Polar Magnitude of filter freq resp
Private filterP(NFFT) As Double   'Polar Phase of filter freq resp

Private IOverlap(NFFT - FILTERTAPS - 1) As Double  'Overlap prev FFT/IFFT
Private QOverlap(NFFT - FILTERTAPS - 1) As Double  'Overlap prev FFT/IFFT

Public AGCHang As Long   'AGC AGCHang time factor

Код парсинга входного буфера inBuffer() показан на рис. 3. Левый и правый сигнальные входы должны быть разобраны на I и Q сигнальные каналы, прежде чем они будут представлены на входе FFT. 16-битовые целочисленные выборки левого и правого каналов во входном буфере перемежаются, поэтому код, показанный на рис. 3, должен использоваться для разделения сигналов. Массивы RealIn() и RealOut() используются для хранения векторов I-сигнала, а массивы ImagIn() и ImagOut() используются для хранения векторов Q-сигнала. Это соответствует номенклатуре, используемой в сложном алгоритме FFT. Не критично, какой из I и Q каналов идет на какой вход, потому что можно просто переставить массивы в коде рис. 3, если боковые полосы инвертированы.

Erase RealIn, ImagIn

For S = 0 To CAPTURESIZE - 1 Step 2    'Копирование I в RealIn и Q в ImagIn
   RealIn(S \ 2) = inBuffer(S + 1)
   ImagIn(S \ 2) = inBuffer(S)
Next S

Рис. 3. Парсинг входных буферов для получения векторов сигналов I и Q.

[FFT: преобразование в домен частот (спектр)]

В части 1 [18] этой серии статей обсуждается, как FFT используется для преобразования дискретных по времени сигналов из области времени (time domain) в область частот (frequency domain, или спектр). Реализация FFT "в лоб" довольно сложна математически и несколько утомительна для кодирования. К счастью, Intel предоставила код с оптимизированной производительностью в виде DLL, который можно вызвать из одной строки кода для этой и других важных функций DSP (см. примечание [4]).

FFT фактически состоит из ряда очень узких полосовых фильтров, выходные сигналы которых называются бинами, как показано на рис. 4. Каждый бин имеет значение магнитуды и фазы, представляющие содержание дискретизированного входного сигнала на центральной частоте соответствующего бина. Перекрытие соседних бинов напоминает выход гребенчатого фильтра, как описано в части 1.

Рис. 4. Выходные "бины" FFT.

PC SDR использует 4096-bin FFT. При частоте дискретизации 44100 Гц полоса пропускания каждого элемента дискретизации составляет 10.7666 Гц (44100/4096), а центральная частота каждого элемента дискретизации равна числу элементов дискретизации, умноженному на полосу пропускания. Отметим на рис. 4, что относительно центральной частоты дискретизированного квадратурного сигнала верхняя боковая полоса расположена в бинах 1 .. 2047, а нижняя боковая полоса расположена в бинах 2048 .. 4095. Ячейка 0 содержит несущую, преобразованную в 0 Гц. FFT, выполняемый на аналитическом сигнале I + jQ, позволяет анализировать положительные и отрицательные частоты раздельно.

Я использовал звуковую карту Turtle Beach Santa Cruz, у которой АЧХ гладкая с точностью 3 dB в диапазоне частот от 10 Гц до 20 кГц (замечание: даташит указывает верхнюю частоту среза 120 кГц, что скорее всего типографическая ошибка, так как максимальная частота дискретизации у этой карты 48 кГц). Поскольку мы оцифровываем RF в квадратурной форме, частота дискретизации фактически удваивается (если 44100 Гц умножить на 2 канала, получается 88200 Гц эффективная частота выборок). Это означает, что выходной спектр FFT будет в два раза больше, чем у одного оцифрованного канала. В нашем случае общая выходная полоса частот FFT составит 10.7666 Гц умножить на 4096 или 44100 Гц. Поскольку большинство звуковых карт дают спад АЧХ 20 кГц, мы вероятно, ограничены общей полосой пропускания примерно 40 кГц.

Рис. 5 показывает вызовы DLL библиотеки Intel для FFT и последующего преобразования векторов сигнала из картезианской (декартовой) системы координат в полярную систему координат. Подпрограмма nspzrFftNip принимает на входе вектора домена времени RealIn() и ImagIn(), и преобразует их в вектора домена частот RealOut() и ImagOut(). Порядок (order) FFT вычисляется в подпрограмме, которая формирует коэффициенты фильтра, что будет обсуждаться далее. NSP_Forw это константа, которая задает для подпрограммы выполнить прямое быстрое преобразование Фурье (FFT).

'FFT:
nspzrFftNip RealIn, ImagIn, RealOut, ImagOut, order, NSP_Forw
'Преобразование Cartesian -> polar:
nspdbrCartToPolar RealOut, ImagOut, M, P, NFFT

Рис. 5. Преобразование сигнала time domain в сигнал frequency domain с помощью FFT.

В декартовой (Cartesian) системе координат сигнал представлен величинами двух векторов: одного на вещественной плоскости или плоскости x, и одного на воображаемой плоскости или плоскости y. Эти векторы могут быть преобразованы в один вектор с величиной (М) и фазовым углом (Р) в полярной системе координат. В зависимости от конкретного алгоритма DSP, который мы хотим выполнить, одна система координат или другая может быть более эффективной. Я использую полярную систему координат для большей части обработки сигнала в этом примере. Подпрограмма nspdbrCartToPolar преобразует выходные данные FFT в полярный вектор, состоящий из значений магнитуды M() и значений фаз P(). Эта функция одновременно выполняет вычисления по формулам Выражения3 и Выражения4 в части 1 [18] этой серии статей.

[Смещение преобразования основной ПЧ к 0 Гц]

Моя оригинальная программа центрировала частоту несущей RF в бине 0 (0 Гц). С помощью этой реализации можно отображать (и слышать) весь спектр 44 кГц в режиме реального времени. Одна из проблем, встречающихся в приемниках с прямым преобразованием или с нулевым ПЧ, заключается в том, что шум увеличивается по существу около 0 Гц. Это вызвано несколькими механизмами: шумом 1/f в активных компонентах, шумом 60/120 Гц (для США, или 50/100 Гц для России) от линий электропередач переменного тока, микрофонным шумом, вызванным механической вибрацией и фазовым шумом гетеродина. Это может быть проблемой для работы со слабым сигналом, потому что большинство людей настраивают сигналы CW для тона 700-1000 Гц. К счастью, большая часть этого шума исчезает выше 1 кГц.

Учитывая, что у нас есть 44 кГц спектра для работы, мы можем сдвинуть цифровую ПЧ на любую частоту в пределах выходного диапазона FFT. Это просто вопрос принятия решения о том, какой бин БПФ обозначить в качестве несущей частоты, а затем сместить гетеродин на соответствующую величину. Затем выполняется копирование соответствующих ячеек для требуемой боковой полосы так, чтобы они располагались на 0 Гц для последующей обработки. В PC SDR я выбрал использовать смещение ПЧ (IF) 11025 Гц, что составляет одну четвертую частоты дискретизации, как показано на рис. 6.

Рис. 6. Диаграмма смещения ПЧ основного диапазона. Частота гетеродина смещается на 11025 кГц, так что требуемая несущая частота сигнала центрирована со смещением 11025 Гц на выходе FFT. Чтобы сдвинуть сигнал для последующей фильтрации, желательные элементы дискретизации просто копируются для центрирования несущей частоты fc на 0 Гц.

На рис. 7 показан исходный код для смещения ПЧ (IF) на 0 Гц. Частота несущей 11025 Гц сдвигается к бину 0, а верхняя полоса частот смещается на бины от 1 до 1023. Нижняя полоса частот смещается на бины от 3072 до 4094. Код позволяет разрешить или запретить смещение IF, как это требуется для передачи.

IFShift = True                     'Принудительно True для demo

   If IFShift = True Then          'Смещение боковых полос от 11025 кГц ПЧ
      For S = 0 To 1023
         If USB Then
            M(S) = M(S + 1024)     'Перемещение верхней полосы на 0 Гц
            P(S) = P(S + 1024)
         Else
            M(S + 3072) = M(S + 1) 'Перемещение нижней полосы на 0 Гц
            P(S + 3072) = P(S + 1)
         End If
      Next
   End If

Рис. 7. Код для преобразования с понижением частоты со смещенной ПЧ основной полосы до 0 Гц.

[Выбор боковой полосы]

Итак, как мы выбираем боковую полосу? Мы храним нули в бинах, которые не хотим слышать. Насколько это просто? Если бы можно было иметь идеальную аналоговую амплитуду и фазовый баланс на дискретизированных входных сигналах I и Q, у нас было бы бесконечное подавление боковой полосы. Поскольку это невозможно, любой дисбаланс будет проявляться в виде изображения в полосе пропускания приемника. К счастью, эти дисбалансы могут быть исправлены с помощью кода DSP либо во временной области до FFT, либо в частотной области после FFT. Эти методы выходят за рамки этого обсуждения, но я могу рассказать о них в следующей статье. Мой прототип, использующий инструментальные усилители INA103, достигает примерно 40 дБ подавления противоположной боковой полосы без коррекции в программном обеспечении.

Код для обнуления противоположной полосы частот предоставлен на рис. 8. Нижняя полоса частот находится в бинах с большими номерами, а верхняя полоса частот находится в бинах с меньшими номерами. Для экономии времени обработки я только обнулил количество бинов, содержащееся в переменной FFTBins.

[Магия фильтрации FFT Fast-Convolution]

Каждый текст DSP, который я прочитал по модуляции и демодуляции одной боковой полосы, описывает подход оцифровки ПЧ (IF). По этому методу A/D-преобразователь оцифровывает сигнал на такой IF, как 40 кГц. Затем сигнал квадратурно преобразуется с понижением частоты на программном обеспечении в основную полосу частот и фильтруется с использованием фильтров с характеристикой finite impulse response (FIR, по-русски этот термин называется фильтр с конечной импульсной характеристикой КИХ) [9]. Такие системы описал Doug Smith QEX в статье "Signals, Samples, and Stuff: A DSP Tutorial (Part 1)" [10]. При таком подходе вся обработка выполняется во временной области (time domain).

Для PC SDR я решил использовать совершенно другой подход, называемый фильтрацией быстрой свертки FFT, FFT fast-convolution filtering (также называемой FFT convolution, или сверткой FFT), который выполняет все функции фильтрации в частотной области [11]. FIR-фильтр выполняет свертку входного сигнала с импульсной характеристикой фильтра во временной области. Свертка - это математическое средство объединения двух сигналов (например, входного сигнала и импульсной характеристики фильтра) для формирования третьего сигнала (отфильтрованного выходного сигнала) [12]. Подход временной области очень хорошо работает для небольшого количества отводов фильтра. Что делать, если мы хотим построить очень высокопроизводительный фильтр с 1024 или более ответвлениями? Время на обработку данных FIR-фильтра может оказаться чрезмерно большим. Оказывается, что важным свойством преобразования Фурье является то, что свертка во временной области равна умножению в частотной области. Вместо прямой свертки входного сигнала с импульсной характеристикой фильтра с оконной структурой, как в случае FIR-фильтра, мы берем соответствующие FFT входного сигнала и импульсной характеристики фильтра и просто перемножаем их, как показано на рис. 9. Чтобы вернуться к временной области, мы выполняем IFFT продукта. Свертка FFT часто быстрее, чем прямая свертка для ядер фильтров длиной более 64 ответвлений, и дает точно такой же результат.

Рис. 9. Блок-схема фильтрации FFT fast-convolution. Коэффициенты импульсной характеристики фильтра сначала преобразуются в частотную область с использованием FFT и сохраняются для повторного использования подпрограммой фильтра. Каждый сигнальный блок преобразуется с помощью FFT и затем умножается на амплитуды частотной характеристики фильтра. Полученный отфильтрованный сигнал преобразуется обратно во временную область с использованием обратного БПФ (IFFT). Подпрограмма наложения/сложения (Overlap/Add) корректирует сигнал для круговой свертки.

Для меня свертку FFT (FFT convolution) легче понять, чем прямую свертку (direct convolution), потому что я мысленно визуализирую фильтры в частотной области. Как описано в части 1 [18] этой серии статей, выходной сигнал комплексного FFT можно рассматривать как длинный набор узких полосовых фильтров, выровненных вокруг несущей частоты (бин 0), что показано на рис. 4. Рис. 10 иллюстрирует процесс свертки FFT преобразованной импульсной характеристики фильтра с преобразованным входным сигналом. Как только сигнал преобразуется обратно во временную область через IFFT, мы должны выполнить процесс, называемый методом перекрытия/добавления (overlap/add). Это происходит потому, что процесс свертки производит выходной сигнал, который равен по длине сумме входных выборок плюс отводов фильтра минус один. Я не буду пытаться объяснить концепцию здесь, потому что она лучше всего описана в ссылках [13].

Рис. 10. Выход фильтрации FFT fast convolution. Когда коэффициенты filter-magnitude умножаются на значения signal-бин, результирующие выходные бины содержат значения только в пределах полосового фильтра.

Рис. 11 предоставляет исходный код для формирования коэффициентов полосового фильтра частотной области. Подпрограмма CalcFilter принимает в качестве параметров нижнюю частоту среза FLow и верхнюю частоту среза FHigh, чтобы сформировать АЧХ фильтра. Затем частоты среза преобразуются в их соответствующие доли частоты дискретизации для использования подпрограммой формирования фильтра nspdFirBandpass. Порядок FFT также определяется этой подпрограммой, основываясь на размере FFT, NFFT. Подпрограмма nspdFirBandpass вычисляет импульсную характеристику полосового фильтра от Fl() до Fh() и длины FILTERTAPS. Затем результат помещается в переменную массива Rh(). NSP_WinBlackmanOpt вызывает формирование окна импульсной характеристики с помощью оконной функции Blackman. Что такое оконная функция, см. DSP Guide [14]. Значение 1, которое передается подпрограмме, приводит к нормализации результата.

Public Static Sub CalcFilter(FLow As Long, FHigh As Long)
Static Rh(NFFT) As Double   'Impulse response for bandpass filter
Static Ih(NFFT) As Double   'Imaginary set to zero
Static reH(NFFT) As Double  'Real part of filter response
Static imH(NFFT) As Double  'Imaginary part of filter response

Erase Ih

Fh = FHigh / Fs             'Вычисление верхней и нижней частот среза
Fl = FLow / Fs              ' как доли частоты дискретизации Fs
BinSize = Fs / NFFT         'Вычисление размера бина FFT в Гц

FFTBins = (FHigh / BinSize) + 50   'Количество бинов FFT в полосе фильтра

order = NFFT                'Вычисление порядка NFFT как степени 2
Dim O As Long

For O = 1 To 16             'Вычисление порядка фильтра
   order = order \ 2
   If order = 1 Then
      order = O
      Exit For
   End If
Next

'Вычисление коэффициентов полосового IIR-фильтра
' с оконной функцией Blackman:
nspdFirBandpass Fl, Fh, Rh, FILTERTAPS, NSP_WinBlackmanOpt, 1

'Вычисление комплексных значений полосового фильтра в частотной области:
nspzrFftNip Rh, Ih, reH, imH, order, NSP_Forw
nspdbrCartToPolar reH, imH, filterM, filterP, NFFT

End Sub

Рис. 11. Код для генерации коэффициентов полосового фильтра в частотной области (frequency domain).

Затем импульсная характеристика преобразуется в частотную область с помощью nspzrFftNip. Входными параметрами являются Rh() вещественная часть импульсной характеристики, и Ih() мнимая часть, которая была установлена на ноль. NSP_Forw предписывает подпрограмме выполнить прямое FFT. Затем мы преобразуем результат частотной области FFT, reH() и imH() в полярную форму с помощью подпрограммы nspdbrCartToPolar. Массивы магнитуд фильтра filterM() и фаз фильтра filterP() сохраняются для использования в фильтре быстрой свертки FFT. За исключением случаев, когда мы вручную изменяем выбор полосового фильтра, отклик фильтра не изменяется. Это означает, что мы должны рассчитать отклик фильтра только один раз, когда фильтр впервые выбран пользователем.

На рис. 12 представлен код для фильтра быстрой свертки БПФ (FFT fast-convolution filter). Используя подпрограмму nspdbMpy2, двоичные коды амплитуды спектра сигнала M() умножаются на двоичные коды амплитуды частотного отклика фильтра filterM() для генерирования результирующих отфильтрованных по месту двоичных кодов амплитуды отклика M(). Затем мы используем nspdbAdd2 для добавления двоичных кодов фазы сигнала P() к двоичным кодам фазы фильтра filterP(), при этом результат сохраняется по месту в отфильтрованных двоичных кодах фазы P().

nspdbMpy2 filterM, M, NFFT   'Умножение бинов магнитуды
nspdbAdd2 filterP, P, NFFT   'Добавление бинов фазы

Рис. 12. Код фильтрации быстрой свертки БПФ (FFT fast convolution filtering) с использованием полярных векторов.

Обратите внимание, что свертка БПФ также может быть выполнена в декартовых (cartesian) координатах с использованием метода, показанного на рис. 13, хотя этот метод требует больше вычислительных ресурсов. Другие применения значений амплитуды частотной области включают в себя усреднение FFT, цифровое подавление (digital squelch) и отображение спектра.

'Вычисление: RealIn(s) = (RealOut(s) * reH(s)) - (ImagOut(s) * imH(s))
nspdbMpy3 RealOut, reH, RealOut_1, NFFT
nspdbMpy3 ImagOut, imH, ImagOut_1, NFFT
nspdbSub3 RealOut_1, ImagOut_1, RealIn, NFFT 'RealIn for IFFT

'Вычисление: ImagIn(s) = (RealOut(s) * imH(s)) + (ImagOut(s) * reH(s))
nspdbMpy3 RealOut, imH, RealOut_2, NFFT
nspdbMpy3 ImagOut, reH, ImagOut_2, NFFT
nspdbAdd3 RealOut_2, ImagOut_2, ImagIn, NFFT 'ImagIn for IFFT

Рис. 13. Альтернативный код фильтрации быстрой свертки БПФ с использованием декартовых (cartesian) векторов.

На рис. 14 показан действительный спектральный выходной сигнал 500-Гц фильтра, использующего широкополосный входной шум и усреднение FFT сигнала в течение нескольких секунд. Это обеспечивает хорошую картину частотной характеристики и формы фильтра. Коэффициент формы 2048-отводного фильтра составляет 1.05 от 3-dB до 60-dB точек (большинство производителей измеряют от 6 дБ до 60 дБ, это более мягкая спецификация). Обратите внимание, что затухание в полосе частот подавления превышает 120 дБ примерно на 250 Гц от точек 3-dB. Это очень высококачественный фильтр!

Рис. 14. Вид АЧХ фактического CW-фильтра на 500-Гц. Фильтрация быстрой свертки БПФ (FFT fast-convolution filtering) используется с 2048 отводами фильтра для получения коэффициента формы характеристики 1.05 от 3 dB до 60 dB вниз и более 120 dB затухания в полосе подавления частот всего 250 Гц за пределами точек 3 dB.

Интересным фактом этого способа является то, что окно применяется к импульсной характеристике фильтра, а не к входному сигналу. Отклик фильтра нормализуется, так что сигналы в полосе пропускания не ослабляются в частотной области. Я полагаю, что эта нормализация отклика фильтра устраняет обычное ослабление, связанное с оконной обработкой сигнала перед выполнением FFT. Чтобы преодолеть такое затухание окна, типично применять перекрытие 50-75% в процессе выборки во временной области и усреднять FFT в частотной области. Буду признателен за комментарии знающих читателей по этой гипотезе.

[IFFT и Overlap/Add - преобразование обратно в Time Domain]

Перед возвратом во временную область (time domain) мы должны сначала преобразовать данные обратно в декартовы (Cartesian) координаты с помощью подпрограммы nspdbrPolarToCart, как показано на рис. 15. Затем путем установки флага NSP_Inv выполняется обратная БПФ (IFFT) подпрограммой nspzrFftNip, которая помещает time-domain выходные данные в массивы RealOut() и ImagOut(). Как уже обсуждалось ранее, мы должны снова реализовать перекрытие с добавлением (overlap and add) порции сигнала из предыдущего цикла захвата, как это описывается в DSP Guide (см. примечание [13]). Ioverlap() и Qoverlap() сохраняют синфазный (in-phase) и квадратурный (quadrature) сигналы с перекрытием от последнего прохода для добавления к новому блоку сигнала с помощью подпрограммы nspdbAdd3.

'Преобразование из полярных значений в декартовы (polar to cartesian)
   nspdbrPolarToCart M, P, RealIn, ImagIn, NFFT

   'IFFT для обратного преобразование во временную область (time domain)
   nspzrFftNip RealIn, ImagIn, RealOut, ImagOut, order, NSP_Inv

   'Overlap и Add от последнего FFT/IFFT: RealOut(s) = RealOut(s) + Overlap(s)
   nspdbAdd3 RealOut, IOverlap, RealOut, FILTERTAPS - 2
   nspdbAdd3 ImagOut, QOverlap, ImagOut, FILTERTAPS - 2

   'Сохранение Overlap для следующего прохода
   For S = BLKSIZE To NFFT - 1
      IOverlap(S - BLKSIZE) = RealOut(S)
      QOverlap(S - BLKSIZE) = ImagOut(S)
   Next

Рис. 15. Код обратного БПФ (Inverse FFT, IFFT) и перекрытия/добавления (overlap/add).

[Цифровая АРУ с переменным временем удержания]

Код цифровой АРУ (AGC) на рис. 16 предоставляет управление усилением с быстрой атакой (fast-attack, быстрое нарастание усиления) и задержкой затухания усиления (decay) и переменным временем удержания усиления (Variable Hang Time). Как атака усиления, так и затухание усиления обе происходят за приблизительно 1 мс, однако время удержания затухания усиления (hang time) может быть установлено в любое желаемое значение с шагом 46 мс. Я решил реализовать атаку/затухание с линейной функцией нарастания, вместо экспоненциальной функции, как описано в коммуникационных текстах DSP [15]. Она работает очень хорошо, и интуитивно понятна для программирования. Блок-схема на рис. 17 описывает логику, используемую в алгоритме АРУ.

If AGC = True Then
   'Если true, то инкремент счетчика AGCLoop, иначе сброс его в 0
   AGCLoop = IIf(AGCLoop < AGCHang - 1, AGCLoop + 1, 0)
   nspdbrCartToPolar RealOut, ImagOut, M, P, BLKSIZE 'Огибающая полярной магнитуды
                                                     ' (Envelope Polar Magnitude)
   Vpk = nspdMax(M, BLKSIZE)                         'Получение пика магнитуды

   If Vpk < > 0 Then                       'Проверка для устранения деления на 0
      G(AGCLoop) = 16384 / Vpk             'Коэффициент усиления AGC с запасом 6 dB
      Gain = nspdMin(G, AGCHang)           'Поиск снижения пикового усиления
   End If

   If Gain > MaxGain Then Gain = MaxGain   'Ограничение усиления до MaxGain

   If Gain < PrevGain Then                 'Усиление AGC уменьшается
      GainStep = (PrevGain - Gain) / 44    '44 Sample ramp = 1 ms attack time
      For S = 0 To 43                      'Скат усиления с периодом 1 мс
         M(S) = M(S) * (PrevGain - ((S + 1) * GainStep))
      Next
      For S = 44 To BLKSIZE - 1            'Умножение остальной огибающей на усиление
         M(S) = M(S) * Gain
      Next
   Else
      If Gain > PrevGain Then              'Усиление AGC увеличивается
         GainStep = (Gain - PrevGain) / 44 '44 Sample ramp = 1 ms decay time
         For S = 0 To 43                   'Нарастание усиления с периодом 1 мс
            M(S) = M(S) * (PrevGain + ((S + 1) * GainStep))
         Next
         For S = 44 To BLKSIZE - 1         'Умножение остальной огибающей на усиление
            M(S) = M(S) * Gain
         Next
      Else
         nspdbMpy1 Gain, M, BLKSIZE        'Умножение огибающей на усиление AGC
      End If
   End If

   PrevGain = Gain                         'Сохранение усиления для следующего цикла
   nspdbThresh1 M, BLKSIZE, 32760, NSP_GT  'Жесткий ограничитель для предотвращения
                                           ' переполнения
End If

Рис. 16. Код цифровой АРУ (Digital AGC).

Рис. 17. Блок-схема алгоритма цифровой АРУ (Digital AGC).

Для последующего описания давайте рассмотрим код на рис. 16 и его алгоритм на рис. 17. Сначала мы проверяем, что АРУ включена. Если это так, то мы инкрементируем AGCLoop, счетчик циклов времени задержки АРУ (AGC hang-time loops). Каждый проход по этому коду соответствует времени задержки 46 мс. PC SDR предоставляет настройки времени hang-time 3 (быстрое АРУ, восстановление 132 мс), 5 (среднее время восстановления, 230 мс), 7 (медленное время восстановления, 322 мс) и 22 (долгое время восстановления, 1.01 секунды). Настройка hang-time сохранена в переменной AGCHang. Как только счетчик hang-time сбрасывается, затухание происходит с линейным наклоном 1 мс.

Чтобы определить требование к усилению АРУ, мы должны детектировать огибающую демодулированного сигнала. Это просто осуществить преобразованием координат из декартовых (Cartesian) в полярные (polar). Значение M() это огибающая, или магнитуда сигнала. Вектор фазы в контексте АРУ может быть проигнорирован. Однако нам потребуется сохранить значения фаз для последующего преобразования в декартовы координаты. Как только у нас есть величины, хранящиеся в M(), становится простым нахождение пиковой величины и сохранение ее в Vpk с функцией nspdMax. После проверки, чтобы предотвратить ошибку деления на ноль, мы вычисляем коэффициент усиления относительно 50% от значения полной шкалы. Это обеспечивает запас в 6 дБ от пика сигнала до полной шкалы выходного значения DAC (ЦАП). При каждом проходе коэффициент усиления сохраняется в массиве G(), чтобы мы могли найти снижение пикового усиления в течение периода зависания с помощью функции nspdMin. Затем коэффициент уменьшения пикового усиления сохраняется в переменной Gain. Обратите внимание, что Gain сохраняется как коэффициент, а не в децибелах, так что преобразование log/antilog не требуется.

Следующим шагом является ограничение Gain значением MaxGain, которое может быть установлено пользователем. Эта система функционирует так же, как управление усилением IF, позволяющее Gain изменяться от отрицательных значений до настройки MaxGain. Хотя такое не предусмотрено в примере кода, это простая задача, чтобы создать элемент управления передней панели в Visual Basic, чтобы вручную установить значение MaxGain.

Далее определяем, нужно ли увеличивать, уменьшать или оставлять без изменений усиление. Если Gain меньше, чем PrevGain (это настройка Gain из блока сигналов, сохраненного на последнем проходе через код), мы линейно уменьшаем усиление по 44 выборкам. Это дает время атаки приблизительно 1 мс при частоте дискретизации 44100 Гц. GainStep - это наклон линейного изменения за время выборки, рассчитанный по значениям PrevGain и Gain. Затем мы постепенно уменьшаем первые 44 выборки на значение GainStep. После перехода к новому значению Gain оставшиеся выборки умножаются на фиксированное значение Gain.

Если Gain увеличивается от значения PrevGain, процесс просто реверсируется. Если коэффициент усиления не изменился, все выборки умножаются на текущую настройку коэффициента усиления. После обработки сигнального блока Gain сохраняется в PrevGain для следующего сигнального блока. Наконец, nspdbThresh1 реализует жесткий ограничитель примерно на максимальном выходном уровне DAC, чтобы предотвратить переполнение выходных буферов с целочисленной переменной.

[Отправка демодулированного или модулированного сигнала в выходной буфер]

Последним шагом является форматирование обработанного сигнала для вывода в ЦАП. При приеме сигнал RealOut() копируется выборка за выборкой, как в левый, так и в правый каналы. Для передачи RealOut() копируется в правый канал, а ImagOut () копируется в левый канал DAC. Если желателен бинауральный прием, сигнал I и Q может быть необязательно послан в правый и левый каналы соответственно, так же как в режиме передачи.

[Управление демонстрационным кодом]

Демонстрационный код SDR (см. примечание [3]) имеет несколько выбранных кнопок для установки времени задержки AGC (АРУ), выбора фильтра и выбора боковой полосы. Код для этих функций показан на рис. 18. Код понятен и легко модифицируется для дополнительных фильтров, разного времени задержки АРУ и других режимов работы. Не стесняйтесь экспериментировать.

Private Sub cmdAGC_Click(Index As Integer)
   MaxGain = 1000   'Maximum digital gain = 60dB
   Select Case Index
   Case 0
      AGC = True
      AGCHang = 3           '3 x 0.04644 sec = 139 ms
   Case 1
      AGC = True
      AGCHang = 7           '7 x 0.04644 sec = 325 ms
   Case 2
      AGC = False           'AGC Off
   End Select
End Sub

Private Sub cmdFilter_Click(Index As Integer)
   Select Case Index
   Case 0
      CalcFilter 300, 3000  'Фильтр 2.7 кГц
   Case 1
      CalcFilter 500, 1000  'Фильтр 500 Гц
   Case 2
      CalcFilter 700, 800   'Фильтр 100 Гц
   End Select
End Sub

Private Sub cmdMode_Click(Index As Integer)
   Select Case Index
   Case 0                   'Изменение режима на USB
      SSB = True
      USB = True
   Case 1                   'Изменение режима на LSB
      SSB = True
      USB = False
   End Select
End Sub

Рис. 18. Код управления для демонстрационной передней панели.

[Полнофункциональная программа SDR-1000]

В SDR-1000, моей реализации PC SDR, содержится значительное количество кода, не проиллюстрированного здесь. Я решил сфокусировать эту статью на основном коде DSP, необходимом для модуляции и демодуляции в частотной области. Как позволит время, я надеюсь написать будущие статьи, которые углубятся в другие интересные аспекты дизайна программного обеспечения.

Рис. 19 показывает завершенную переднюю панель SDR-1000. Большинство функций пользовательского интерфейса интуитивно понятны. Вот несколько интересных возможностей SDR-1000:

• Отображение спектра сигнала в реальном времени, с настройкой частоты одним щелчком мыши.
• Сдвоенные независимые VFO с считыванием базы данных о выделении полосы частот. Пользователь может легко получить доступ к базе данных планирования полосы частот и изменить ее.
• Настройка колесом мыши с возможностью изменения скорости настройки щелчком на колесе.
• Многофункциональный цифровой и аналоговый измеритель мгновенного и среднего уровня сигнала, усиления АРУ, входного сигнала АЦП и уровней выходного сигнала ЦАП.
• Обширное управление VFO, диапазоном и режимом. Кнопки переключения диапазонов также обеспечивают многоуровневую память в том же диапазоне. Это означает, что при многократном нажатии на кнопку данной полосы она будет циклически проходить через последние три частоты, посещенные в этом диапазоне.
• Практически неограниченные возможности памяти предоставляются через интерфейс базы данных Microsoft Access. Память включает в себя все ключевые настройки радиостанции по частоте. Частоты также могут быть сгруппированы для сканирования.
• На передней панели предусмотрено десять стандартных настроек фильтра, а также независимые непрерывно изменяемые фильтры как для CW, так и для SSB.
• Локальные и UTC часы реального времени.

Учитывая возможности Visual Basic, возможность расширения пользовательского интерфейса практически безгранична.

Рис. 19. Передняя панель программы SDR-1000.

В коде DSP можно сделать гораздо больше для настройки PC SDR данного приложения. Например, Leif Åsbrink SM5BSZ делает под Linux интересную обработку слабого сигнала, отражаемого от поверхности Луны [16].

Кроме того, Bob Larkin W7PUA использует DSP-10, который он впервые описал в сентябрьских, октябрьских и ноябрьских выпусках QST за 1999 год, для экспериментов со слабым сигналом, принимаемым за горизонтом микроволнового распространения [17].

[Что будет в последней статье цикла]

В заключительной статье я планирую описать текущую разработку аппаратного обеспечения SDR-1000. Будет включен анализ компромисса между распределением усиления, показателем шума и динамическим диапазоном. Я также рассмотрю различные подходы к аналоговому AGC и рассмотрю управление частотой с использованием AD9854 квадратурной DDS.

[Ссылки и примечания]

1. G. Youngblood, AC5OG, "A Software Defined Radio for the Masses, Part 1", QEX, Jul/Aug 2002, pp 13-21.
2. G. Youngblood, AC5OG, "A Software Defined Radio for the Masses, Part 2", QEX, Sep/Oct 2002, pp 10-18.
3. Демонстрационный код этого проекта можно загрузить с ARRLWeb по ссылке www.arrl.org/qexfiles/. Ищите файл архива 1102Youngblood.zip.
4. Функции библиотеки Intel Signal Processing Library предоставлены в пакете Intel Performance Primitives (IPP, Version 3.0, beta) для процессоров Pentium и архитектур Itanium. Оценочная копия IPP доступна бесплатно для загрузки по ссылке developer.intel.com/software/products/ipp/ipp30/index.htm. Коммерческое использование IPP требует полной лицензии. Не используйте IPP с демонстрационным кодом, поскольку он был протестирован только на предыдущей библиотеке обработки сигналов.
5. D. Hershberger, W9GR, and Dr S. Reyer, WA9VNJ, “Using The LMS Algorithm For QRM and QRN Reduction,” QEX, Sep 1992, pp 3-8.
6. D. Hall, KF4KL, “Spectral Subtraction for Eliminating Noise from Speech”, QEX, Apr 1996, pp 17-19.
7. J. Bloom, KE3Z, “Correlation of Sampled Signals”, QEX , Feb 1996, pp 24-28.
8. R. Lyons, Understanding Digital Signal Processing (Reading, Massachusetts: Addison-Wesley, 1997) pp 133, 330-340, 429-430.
9. D. Smith, KF6DX, Digital Signal Processing Technology (Newington, Connecticut: ARRL, 2001; ISBN: 0-87259-819-5; Order #8195) pp 4-1 through 4-15.
10. D. Smith, KF6DX, “Signals, Samples and Stuff: A DSP Tutorial (Part 1)”, QEX (Mar/Apr 1998), pp 5-6.
11. Информацию по свертке БПФ (FFT convolution) можно найти по следующим ссылкам: R. Lyons, Understanding Digital Signal Processing, (Addison-Wesley, 1997) pp 435-436; M. Frerking, Digital Signal Processing in Communication Systems (Boston, Massachusetts: Kluwer Academic Publishers) pp 202-209; and S. Smith, The Scientist and Engineer’s Guide to Digital Signal Processing (San Diego, California: California Technical Publishing) pp 311-318.
12. S. Smith, The Scientist and Engineer’s Guide to Digital Signal Processing (California Technical Publishing) pp 107-122. Это доступно для свободной загрузки на сайте www.DSPGuide.com.
13. Overlap/add method: Ibid, Chapter 18, pp 311-318; M. Freirking, pp 202-209.
14. S. Smith, Chapter 9, pp 174-177.
15. M. Frerking, Digital Signal Processing in Communication Systems, (Kluwer Academic Publishers) pp 237, 292-297, 328, 339-342, 348.
16. See Leif Åsbrink’s, SM5BSZ, веб-сайт ham.te.hik.se/homepage/sm5bsz/.
17. See Bob Larkin’s, W7PUA, домашняя страничка www.proaxis.com/~boblark/dsp10.htm.
18. A Software-Defined Radio for the Masses, часть 1.
19. A Software-Defined Radio for the Masses, часть 2.
20. A Software-Defined Radio for the Masses, часть 4.