Спектральное послание, взято из ролика на этой страничке.

Кусок СИБИ и 10-тки. http://radiokoteg.kiev.ua/RS/FREQ/CB/2015/HDSDR_20151101_124830Z_27771kHz_RF.wav

! 1.9 ГИГ!

Народ, подскажите, как правильно вывести IQ в виде спектра? Я что-то запутался.

Есть HackRF, с него читаю IQ отсчеты. Запускаю RX, получаю callback в которых лежит byte-массив на 260Кб. Это последовательно IQ-отсчеты по 8 бит каждый. Что мне сделать чтобы вывести на экран в виде спектра частот?

Я рассматриваю отсчеты как последовательно I/Q, применяю к ним комплексное FFT. Далее беру абсолютные значения из комплексного результата, вывожу на экран. Спектр вижу, но он почему-то развернут наоборот, и местами есть лишние станции (сравнивал с gnuradio).
Чувствую что я что-то упускаю, но не пойму где.

DVE
Спектр можно развернуть наоборот, если поменять местами все I-отсчёты с Q-отсчётами (переставить массивы I и Q перед FFT).

А с происхождением"зеркалок", может оказаться сложнее разобраться. Не исключено, что они уже есть в исходных IQ-данных из-за какого-то разбаланса между I- и Q-каналом; в профессиональной проге может быть есть какой-то алоритм коррекции IQ-разбаланса, а в самодельной проге нет и "зеркалки" видны.

( P.S. Если скинете сюда ваш массив (желательно в формате wav со стандартным заголовком; или в текстовом файле .txt с двумя столбиками: столбик I-отсчётов и столбик Q-отсчётов в виде десятичных чисел, причём каждая пара чисел в строке пусть разделяется пробелом), то смогу показать Вам спектр и спектрограмму ваших данных в разных режимах обработки в своей самоделке для сравнения с вашим результатом, и если будет надо, то и обсудим алгоритм. Аналогично можно на совсем простых тестовых массивах сделать сравнительную проверку вашего алгоритма изображения спектра.

Другой вариант: посмотрите спектрограммы ваших IQ-данных в проге SA или HDSDR; если будут видны те же "зеркалки", то придётся гуглить статьи о коррекции IQ-баланса в исходных данных, или вникать в исходники gnuradio )

Спасибо, вечером скину, интересно сравнить.

Алгоритма пока особого нет: я беру массив пар чисел I и Q, применяю к ним FFT, получаю комплексный массив. Далее тупо вывожу на экран вычисляя значения log10(i*i + q*q). Спектр при этом кстати перевернут.
Я что-то забываю, или в теории этого достаточно?

Когда я смотрю спектр в gnuradio, я вижу его в обе стороны от нуля, в комплексной плоскости. Я вычисляю абсолютные значения, может из-за этого левая и правая часть накладываются в виде "зеркалок"?

Мысль про SA хорошая, но у меня Mac OS :)

Переворот спектра это ерунда; чистая условность: что в двухканальном сигнале считать I-отсчётами, а что Q-отсчётами; подавайте на вход комплексного FFT массивы Q и I вместо I и Q, и спектр перевернётся.

Я что-то забываю, или в теории этого достаточно?
Вообще-то, не достаточно. Надо до FFT отсчёты I[n] и Q[n] в каждом канале умножить на коэффициенты W[n], называемые "окном", которые постепенно уменьшаются до нуля (или почти до нуля) к началу и к концу массива. Это так, если массивы не очень большие, и FFT применяется сразу ко всем входным I[n] и Q[n].

Но лучше исходные массивы I[n] и Q[n] разбить на несколько блоков одинаковой длины. В каждом блоке числа I[n] и Q[n] умножить на окно W[n], как если бы блок был самостоятельным сигналом. Затем к каждому блоку применить FFT. И вычислить в каждом блоке "спектр мощности", т.е. массив

P[k]= Re[k]*Re[k]+Im[k]*Im[k]

для каждого комплексного числа на выходе FFT; k - номер комплексного числа в выходном массиве FFT (называется "номер бина"), Re и Im - означают операции взятия действительной и мнимой части комплексного числа.

Дальше всё зависит от того, хотите ли Вы изобразить график спектра всего сигнала, или хотите спектрограмму. Если хотите получить график спектра, то все значения P[k] от каждого блока при одном и том же k надо просуммировать (это называется "некогерентное сложение спектров" блоков, и теория учит, что при этом немножко улучшается отношение "сигнал/шум" на итоговом графике: пики сигналов будут более различимы на фоне "травы").

Брать или не брать логарифм от получившихся суммарных значений P[k] - дело вкуса: с логарифмом кратинка будет менее контрастная, шумы и пики слабеньких сигналов (а также боковые части спектра АМ-сигналов) на ней будут видны примерно наравне с пиками сильных сигналов (и с пиками несущих). Без логарифма картина более контрастная; как промежуточный вариант (мне он нравится больше всего, просто по привычке) - можно извлечь квадратный корень из P[k].

А если хотите получить спектрограмму, т.е. картину типа "водопад", то не надо суммировать P[k] от блоков, а надо изобразить для каждого номера k число P[k] (или его логарифм, или корень из него, - по вкусу) "цветной точкой". Номер k служит номером частоты на шкале частот спектрограммы, а номер блока служит номером момента времени на оси времени спектрограммы.

Да, забыл сказать самое главное: блоки следует выбирать перекрывающимися; наименьшее приемлемое перекрытие - половинное: начальный номер данных следующего блока находится в середине предыдущего. Можно и большее перекрытие делать. Т.е. в проге желательна возможность пользователю задавать руками длину блока FFT и длину шага (что равноценно заданию перекрытия блоков) - смещение начального номера данных от блока к блоку.

Если пользователь выбирает очень длинные блоки, то на спектрограмме улучшается "разрешение по частоте", но замываются детали изменения сигналов во времени. Если блоки короткие, то - наоборот. При удачном выборе длины блока и шага спектрограмма не сильно изменяется при изменении этих параметров вблизи выбранных значений - и это есть признак того, что в неё можно верить. Если же картинка сильно меняется при изменении параметров обработки данных, значит результат характеризует в большей мере алгоритм обработки, нежели информацию, содержащую в исходных данных. Примерно так как-то...

А ещё вместо "водопада" можно изображать "живой спектр": шкалы времени на нём нет, есть только шкала частот (привязанная к номерам бинов), и на экран выводится график P[k] (или логарифм и т.п.) сначала для одного блока, потом с какой-то маленькой задержкой во времени - график P[k] для второго блока, потом - для третьего и т.д.

блоки следует выбирать перекрывающимися;
С окном Ханна, которое специально предназначено для половинного перекрытия, в таком случае, спектр получается наименее искаженным и может быть преобразован обратно в сигнал с минимальными потерями информации.

Sinus
Спасибо, реально полезно.

Еще тупой вопрос - а как в gnuradio спектр центрирован относительно нуля получается? Мы применяя FFT и Abs, получим на экране только правую половину?

DVE
применяя FFT и Abs, получим на экране только правую половину?
Не обязательно, зависит от реализации модуля FFT. FFT всегда дает на выходе две половины, просто для упрощения процедуры их могут либо сложить, либо одну выбросить, в итоге выходной массив сокращается в 2 раза. Чисто для визуального восприятия, это не принципиально, а для дальнейшей обработки данных потребуются обе половины.

Вообще, FFT игрушка серьезная, вот недавно приобрел "учебное пособие", для дальнейшего изучения возможностей функции на практике. Уже прикручено к стене. :)

DVE
как в gnuradio спектр центрирован относительно нуля получается? Мы применяя FFT и Abs, получим на экране только правую половину?

Что делается в gnuradio, я без понятия. А принцип центровки спектра не очень давно обсуждался в грубых чертах в другой теме; там у меня есть уже нарисованная картинка, а до неё, кстати, там есть маленькие примерчики IQ-сигналов в wav-файлах и картинки спектров и спектрограмм к ним (ссылаюсь на те посты не для саморекламы,- я ведь так и остался дубом в ЦОС и в программировании, - а просто чтобы то немногое, что вроде знаю, не печатать здесь повторно :)

Функция Abs сама по себе никак не влияет на спектр частот, получаемый из входных IQ-данных с помощью FFT; об этом Вы даже и не беспокойтесь: ни Abs, ни log, ни любые другие подобные действия c выходными числами из FFT не испортят спектр.

Проблема может быть в самом "FFT", если Вы пользуетесь непроверенным вами алгоритмом "FFT", ибо бывают разные алгоритмы.

Например, в простеньком пакете Mathcad есть два разных "FFT": первый из них вызывается как функция FFT() - она принимает на входе только один массив действительных чисел, скажем R[n], а на выходе возвращает массив комплексных чисел F[k], из которого затем операциями Re() и Im() мы можем выделить два вдвое более коротких массива действительных чисел. Из них строим массив P[k], и он да, изображает только "правую половину" от нуля на шкале частот. Но "левая половина" здесь и не нужна, потому что теория учит: у действительного массива R[n] "левая половина" спектра мощности (т.е. в области отрицательных частот) зеркально симметрична к "правой половине".

Однако для FFT-обработки IQ-сигнала (который двухканальный!) нам нужна функция, принимающая на входе поотдельности два массива. В Mathcade такая функция называется CFFT(), она принимает на входе массив комплексных чисел I[n]+j*Q[n], что равноценно двум массивам действительных чисел: I[n] и Q[n]. На выходе нужная нам функция должна возвращать тоже два массива действительных чисел: ReF[k] и ImF[k], - вещественную и мнимую части комплексной фурье-амплитуды F[k]. Здесь длина обоих массивов на выходе уже равна длине массивов на входе. Номера k=0, 1, 2, ... нумеруют частоты от "левого конца" спектра до "правого конца" (равного половине частоты дискретизации в IQ-сигнале); равная нулю фурье-частота находится в центре этого интервала.

Теория учит нас: раз входной сигнал был 2-канальным, т.е. это был IQ-сигнал, то его спектр P[k] в области отрицательных фурье-частот не обязан быть зеркальным отражением спектра в области положительных фурье-частот, как было бы для 1-канального ("моно") входного сигнала.

( В Mathcade функция CFFT() возвращает один массив, но зато комплексный, F[k]= ReF[k] +j*ImF[k]. Как организован ввод/вывод FFT в gnuradio или в других sdr-прогах, я не знаю; знаю только, что принцип раскладки данных по массивам должен быть именно такой; но вдруг там программисты сделали ввод/вывод в виде чередующихся Re- и Im- частей, или ещё как-нибудь. Т.е. в принципе не исключён случай, что надо детально копаться в исходнике, а не брать вслепую оттуда готовый код FFT... или уж по крайней мере надо заранее протестировать его. )

DVE
Сейчас я освежил в своей памяти, как организованы данные на выходе CFFT() в Mathcad-e и вспомнил, что она выдаёт номера выходных чисел не такие, какие хочется иметь для шкалы частот в sdr-проге: она сначала (при k=0, 1, 2, ...) выдает спектр положительных фурье-частот (от нулевой частоты при k=0, до Fdiscr/2 при k=N/2, где N - количество отсчётов в I или Q массиве во входном блоке), а дальнейшим номерам k соответствуют отрицательные частоты: от почти -Fdiscr/2 до нулевой частоты при k=N.

Если у Вас алгоритм FFT тоже так выдаёт результат, то наверняка в этом и вся загвоздка. И если так, то тогда перед выводом на экран достаточно перенумеровать массив спектра мощности P[k].

Попробуйте вот такой цикл перенумеровки (хотя лучше бы сначала глянуть на числа в каком-нибудь маленьком тестовом массиве, чтобы понять способ нумерации выходных данных в вашей конкретной реализации FFT):

Пусть N - количество отсчётов в одном канале (I или Q), оно чётное, нумерация элементов массива ведётся от нуля до N-1. Цикл:

Для k от 0 до (N/2)-1 с шагом 1 приравнять:
RenP[k] := P[k+(N/2)]
RenP[k+(N/2)] := P[k]
Конец цикла.

Здесь RenP[k] - это перенумерованный массив, для вывода на экран. В нём номера принимают значения k = 0, 1, 2, ... , N-1. Теория учит, что спектр периодичен, поэтому можно при желании периодически продолжить этот массив, т.е. положить, что: RenP[N] = RenP[0], RenP[N+1] = RenP[1], и так далее.

Сомневаюсь, что это будет кому-то интересно, но чем черт не шутит.) По результатам тестирования "учебного пособия", обозначенного двумя постами выше: Записываем одну струну (для примера, нота "До" малой октавы) и выделяем первые 5 гармоник, с помощью затертого до дыр FFT.



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

Sinus, Programmist
Огромное спасибо за подсказки, поигрался с размером FFT, вроде заработало.

Но периодически на спектре заметны "фантомные" станции, которые при скроллинге дальше исчезают, частота "фантомов" получается = частоте основной станции +/- sample_rate.

PS: Кстати, если кому интересно покрутить IQ-файл от HackRF (записан штатной утилитой):
http://www.filedropper.com/recordiq (30 секунд записи, центр 100МГц).

DVE
Периодические повторы в спектре с периодом sample_rate как раз должны быть, всё правильно.

На шкале частот sdr-приёмника после переработки IQ-данных в картинку спектра смысл имеет только полоса частот шириной = sample_rate. А вне этой полосы идут повторы спектра, их на экран выводить не нужно.

> На шкале частот sdr-приёмника после переработки IQ-данных в картинку спектра смысл имеет только полоса частот шириной = sample_rate.

Да, это понятно, я имел в виду то, что на этой полосе иногда видны станции, которые реально на sample_rate выше.

Это странно. Если такие "фантомы" видны и в других прогах, то, значит, приёмник выдаёт глючные IQ-данные. А если только в вашей самодельной, то где-то есть всё-таки какие-то ошибки у Вас в алгоритме построения спектра, и придётся тест-сигналы погонять - с небольшим количеством маркерных палок в спектре, посмотреть, возникнут ли какие-нибудь глюки.

если кому интересно покрутить IQ-файл от HackRF (записан штатной утилитой):
Жаль, что он не имеет wav-заголовка... :( Такой файл не могу посмотреть, придётся прикручивать вручную header в hex-редакторе, но для этого надо знать sample_rate, разрядность и количество отсчётов на канал; это долго...

Sinus
Не знаю, то ли особенность железа, то ли алгоритмов, то ли так и должно быть.
Вот здесь на картинке программы GQRX видно:

Увеличить

Если взять сильный сигнал и сдвинуться по спектру влево, то "кусок" сигнала проявляется с другой стороны на краю. Если сдвинуться чуть дальше, эффект пропадает. Возможно это как-то давится программно, не разбирался.

DVE
Назовите, пожалуйста, sample_rate в том файле, который Вы залили.

Вроде, я вижу в нём несколько FM-станций и смогу их демодулировать, но мне надо знать значение sample_rate. Если буду знать, то выложу картинку спектра в вашем том файле и тогда обсудим конкретно результаты (скорее всего завтра). Пока могу сказать только, что сигнал катастрофически слабый в файле, и его большая часть - постоянка.

У меня картинка спектра другая, не такая как показанная Вами; у меня чётко видно несколько станций с частотной модуляцией.

8bit, 2000000 sample rate.
Я перезалил уже другой в wav: http://www.filedropper.com/hackrf20160310001939z100000khziqwav

Сильно заморачиваться не надо, идея-то вроде понятна, а с мануальной компенсацией I/Q я заморачиваться не буду пока что. Файл можно посмотреть разве что как образец iq c HackRF.

> Пока могу сказать только, что сигнал катострофически слабый в файле, и его большая часть - постоянка.

Наверно я gain не выставил в настройках записи, в консольной программе-то не видно визуально что пишется. Образец плохой получился :)

> У меня картинка спектра другая, не такая как показанная Вами

Картинка просто как пример того эффекта, когда станция "с другой стороны" спектра видна.

Сейчас картинку для первого сделаю и покажу (а демодуляцию и обсуждение - завтра)

Вот спектрограмма первых 20 секунд 0.5 сек из первого вашего файла. По горизонтали время, по вертикали частота. Центральную частоту (Center Freq) выставил 100 000 000 Hz, как Вы и сказали. Спектр может быть и перевёрнут.

Увеличить

Пять станций тут вижу; а в верхней части рисунка ещё что-то слабое непонятное.

UPD: ошибся я насчёт длительности фрагмента - при дискрете 2 МГц это кусок сигнала в 0.5 сек всего.

Второй файл тоже скачал, но посмотрю завтра вечером (сейчас надо спать идти, завтра на работу вставать... Спокойной ночи :-))

Sinus
Спасибо. С моими результатами вроде совпадает.
Самая нижняя станция на Вашем скрине - 100.9MHz как раз "фантом" которого реально нет.

GQRX этот диапазон показывает так (при бОльшей sample rate):

Увеличить
Как можно видеть, на 100.9 в реале пусто.

DVE
На Вашей записи отчетливо читаются 5 станций, в их числе самая нижняя и самая верхняя, никаких "фантомов" нет, а вот правильность формирования I/Q сигнала вызывает некоторые сомнения.

Если Вы про уровни сигналов, то да - утилита записи консольная, и я неправильно выставил усилениe приемника.

DVE
про уровни сигналов
Нет, я про "центральную частоту", полученную в результате конвертации I/Q в реальный сигнал, имхо, в правильном реальном сигнале никаких центральных частот не должно быть видно.

Да, на HackRF есть пик в центре на нулевой частоте, видимо это аппаратная особенность.

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

Тут не знаю, может быть.
Вы кстати смотрели WAV или первый IQ-файл? Поддержка wav у них экспериментальная и в стадии беты, может что-то и не то.