1. Берём 4096 бина от IQ последовательности.
2. Выдёргиваем несколько бинов (штук 200) из любого подинтервала, куда попадает интересующий нас сигнал.
3. Копируем в новый буфер. Учитывая то, что в буфере сначала должна идти положительная часть, потом начинаться конец "отрицательной".
4. Делаем IFFT, получаем временное представление интересующего нас сигнала, центральная частота которого перенесена "в ноль" (не смог сформулировать куда она там перенесена).

Короче, всё работает. Таким методом можно выдёргивать кучу каналов одновременно после одного FFT исходной IQ последовательности.

Вопрос вот в чём. Если мы попробуем перенести 8 бинов: (4 первых бина в отрицательную часть спектра, 4 в положительную) - всё работает. Если попробуем перенести 6 бинов: 3 в отрицательную, и 3 в положительную, то всё ломается. Всё ломается, если половинка половины становится нечётной. На выходе слышим не могу описать что. Как будто периодически ломается фаза или как будто подмешали периодически идущие дельта-импульсы в сигнале. Условно, на частоте выходного сигнала 48000 через каждые 197 семплов примерно вижу какие-то выбросы по амплитуде в разные стороны)

Хочется понять куда можно покопать, чтобы осознать что происходит.



Картинка. Когда в частях "1" или "2" становится нечётное число бинов (но равное!), то происходит жесть. В целом и в "отрицательную" и в "положительную" часть уходит одинаковое количество энергии, казалось бы.

Происходят фазовые искажения из-за того, что FFT/IFFT вычисляется над непериодическим сигналом. В частном случае с 8 бинами видимо повезло и уровни по краям оказались достаточно низкими, чтобы замаскировать искажения, с 6 бинами не повезло.
Что делать...
Во первых, и FFT , и IFFT нужно делать с перекрытием и умножать на окно как при прямом, так и при обратном преобразовании.
Во вторых, желательно, чтобы число бинов было степенью двойки как в прямом, так и в обратном преобразовании.
В третьих, лучше всего, получается, когда на выходе остается 1 бин, и, соответственно, IFFT делать уже не нужно, как и умножать на окно. Только в этом случае удастся полностью избавиться от искажений.

Да окно (горб синуса) я использую и перекрытие при FFT 50%

А, вы наверное имеете ввиду применение окна при выдирании куска спектра? Сейчас я выдираю кусок спектра, вычисленного от "заоконенного" сигнала, наложенного на свой 50% сдвиг, но получив FFT я выдираю оттуда 200 бинов без их "заоконивания" и вычисляю iFFT, а результат iFFT так же со сдвигами 50% суммирую. Хехе, возможно стоит попробовать оконить результат iFFT, а может оконить 200 бинов...

Так, моя схема.

signal_win = window(signal_src)
bins = fft4096(signal_win)
band = bins[from 3000 to 3200]
bans = window(bins) // ЭТОГО У МЕНЯ НЕТ -- это важно видимо...?
band = make_good(band) // переставить полож и отриц части как надо
signal_output = inverse_fft(band)

pavelkolodin
bins[from 3000 to 3200]
размер тут должен быть степенью двойки

Можно делать так:

1. Сдвигаете окно на 75%
2. Накладываете оконную функцию (для 75% перекрытия с последующим обратным преобразованием подходит окно Ханна) - можете почитать тут: https://habr.com/ru/articles/430536/
3. Переходите в частотный домен (размер 4096)
4. Переставляете как надо две части
5. Вырезаете нужный частотный кусочек размером степенью двойки. Если вам нужен размер 200 - то берете ближайшую степень двойки - т.е. 256, а ненужные 56 бинов просто зануляете.
6. Обратно во временной домен
7. Окно Ханна на полученный кусок
8. Суммируете с перекрытием 75%

> размер тут должен быть степенью двойки

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

pavelkolodin
Ну в общем вы правы. Распространенные реализации БПФ работают для степени двойки.

> Ну в общем вы правы. Распространенные реализации БПФ работают для степени двойки.

На самом деле я делаю вот что.
Мой вход: 1 мегасемпл/сек I/Q 8 битные, 4096 бинов fft.
1000000 / 4096 = 244.140625 раз в секунду я делаю FFT.
Это не учитывая нахлёст, с нахлёстом 50% - в два раза чаще. Остановимся на рассчётах без нахлёста для ясности, с нахлёстом вы просто делается всё это дважды чаще, просто с нахлёстом в 50%. Технология "нахлёст" - это брать входной сигнал по 4096 семпла не впритык, а с шагом 2048 или 1024 (нахлёст 75%), взятое умножать на окно и от этого делать FFT.
Продолжим для варианта "без нахлёста".
244.140625 -- герц внутри одного fft бина.
Сколько бинов занимает 12500-герцовый канал PMR-рации?
12500 / 244.140625 = 51.2 бина, округлим до 52.
Вот эти 52 бина я и выдираю 244.140625 раз в сек -- в них вся информация о нужном канале 12500.
Теперь, во что хочется это преобразовать во временной области?
В сигнал 48000 гц, 1 канал, например.

Итак, со входа входа идут 244.140625 раз в секунду пачки по 4096 бинов (или в 2 раза чаще для нахлёста 50%). Значит надо столько же раз в секунду делать IFFT на какое число бинов, зная ширину бина в герцах?

48000 / 244.140625 = 196.608 бинов должно быть на выходе.
Оставим 196 бинов для простоты. Да, выходной сигнал в итоге получится номинально не 48000, а 47851.5625.
Но сохраним мы его в 48000, может никто и не заметит.
(48000 - 244.140625 * 196) / 2 = ПЛЮС 74.2 герца -- вот такая ошибка в частотах будет на выходе. Все голоса будут звучать выше (тоньше) на 74 гц. Допустим, нам на это пофиг. Не охота возиться с передескретизацией. Можно подобрать другие параметры всей этой комбинации, чтобы минимизировать ошибку при минимуме вычислений (отказываясь от передескретизации).

Короче, мне надо упихать 52 бина в 196 бинов. Как? Наши 52 бина - это же "с нулём посередине". Ноль - это какбэ несущая.То есть, всё что звучит в бине номер (52 / 2) и далее - это частоты выше центра нашего канала 12500, а всё что находится в бине (52 / 2) - 1 и ниже -- это первая половина частот нашего канала 12500 - туда попадает энергия, когда ЧМ-модуляция делает снижения частоты, передавая нижнюю полуволну исходного звукового сигнала.

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

Так вот, получается в наших целевых 196 бинах всё должно лежать так:
[0...52/2] -- сюда кладём вторую половину наших выдернутых 52 бинов.
[196-1...196-(1 + 52/2)] -- сюда кладём первую половину наших 52 бинов.

Делаем IFFT и готово. Можно демодулировать FM-демодулятором. Чтобы выход был без подземных стуков, примени нахлёст с окном.

Так вот мои проблемы начинались, когда половинка этих 52 семплов оказывалась нечётной. Например 52/2 -- это НЕ проблема, потому что 52/2 = 26 - чётное. Если вместо 4096 в этой истории появляются 8192, то жопа отваливается, потому что вместо 52 получается 122, а 122/2 = 61.

То есть, в своей реализации я даже для IFFT не использую степень двойки (использую 196 бинов, например) и "вроде как-то работает", даже шумов особо нет и фаза вроде норм, FM-демодулятор выдаёт ясный приемлимый сигнал.

Сколько бинов занимает 12500-герцовый канал PMR-рации?
Предлагаю взять 1 бин (или лучше 0.5/0.25 бина) уменьшив размер FFT и одним FFT сразу убить 3 зайцев:
1. Нарезать широкую полосу на узкие полосы
2. Понизить частоту дискретизации
3. Отфильтровать полученные полосы

Берем sample rate 1.6 Msps (в допустимых пределах для RTL-SDR), размер FFT 256, перекрытие 75% и получаем полосы по 25 kHz с шагом 6.25 kHz - то, что нужно.

Оставим 196 бинов для простоты. Да, выходной сигнал в итоге получится номинально не 48000, а 47851.5625.
Но сохраним мы его в 48000, может никто и не заметит.
Заметит. Звуковая карта. И очень быстро заметит. Нужно ресэмплировать.

Но FFT хранит выборки так: первая половина - положительные частоты, вторая половина - отрицательные.
Первая половина это поворот фазы от 0 до 180 градусов, вторая половина это поворот фазы от 180 до 360 градусов или от -180 до 0 градусов, что одно и то же. FFT всегда закольцовано и на входе и на выходе. Что тут может быть не понятно.

Делаем IFFT и готово. Можно демодулировать FM-демодулятором.
Если демодулируем 1 канал (ну или 2) то в FFT нет никакого смысла - FFT будет медленнее FIR с прореживанием при сопоставимом качестве.
Если демодулируем много каналов (те же 16 каналов PMR с детьми и строителями), то выбираем размер FFT так, чтобы в один бин умещалась полоса одного канала с запасом раза в 2..3. Тогда можно будет просто брать нужные бины после каждого FFT, дофильтровывать FIR при необходимости с доворотом фазы для убирания постоянного смещения после демодуляции и подавать на FM демодулятор. Так получается и быстрее и проще.
FFT можно делать с перекрытием 50 или 75% - decimation ratio уменьшится соответственно в 2 или 4 раза со значительным улучшением качества фильтрации.
Зачем тут вообще приплетают IFFT я не понимаю.

Фокус с FFT/IFFT имеет смысл в случае когда порядок FIR становится несколько миллионов или хотя бы сотен тысяч и прореживание не требуется. Тогда да, схема "умножение на окно, FFT, умножение на FFT от taps, IFFT, умножение на окно, копирование с перекрытием" начинает выигрывать по скорости у FIR.

P.S. Под поворотом фазы имеется в виду поворот на каждом шаге (инкремент), а не умножение на константу.

> Берем sample rate 1.6 Msps

rtl-sdr такое не умеет, он умеет 250K, 1M, 2M вроде и всё. Ну идею понял, подумаю.

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


> Первая половина это поворот фазы от 0 до 180 градусов, вторая половина это поворот фазы от 180 до 360 градусов или от -180 до 0 градусов, что одно и то же. FFT всегда закольцовано и на входе и на выходе. Что тут может быть не понятно.

С фазой не понял. Ладно, я в принципе слабо понимаю FFT, надо ещё раз вникать. Ведь когда мы кидаем I семплы в real часть комплексного числа, а Q-семплы в Im-часть комплексного числа и делаем над этим FFT и потом получаем частотную картинку от минус половины диапазона частот до плюс половины. Ну, всмысле: "перевёрнутость" синуса (Q-части) несёт в себе информацию о частотах ниже выбранной текущей частоты local oscillator SDR.

Ищу исходник, фурье преобразования. Без использования библиотек. Паскаль, бейсик или си без плюсов и #.

rtl-sdr такое не умеет
Умеет:
https://github.com/osmocom/rtl-sdr/blob/master/src/librtlsdr.c#L1099-L1104

В файл сохраняю, там пофиг.
Тогда можно в заголовок файла, если он есть, записать правильный sample rate.

> Умеет:
> https://github.com/osmocom/rtl-sdr/blob/master/src/librtlsdr.c#L1099-L1104


Не, там же в железе/прошивке выбор между фиксированным набором частот вроде. Оно не может работать на 270000 ГЦ например, только 250000.

http://www.radioscanner.ru/uploader/2023/scrin-2024-10-09-15-.jpg

Не, там же в железе/прошивке выбор между фиксированным набором частот вроде.
Не знаю, откуда скрин, но допустимый диапазон sample rate можно посмотреть в исходниках librtlsdr в указанном месте. Там же можно подсмотреть формулу для вычисления точного значения sample rate. Если у вас в каком-то ПО реализован только выбор из фиксированного набора, то могу только порекомендовать открыть issue в репозитории этого ПО или сразу PR с исправлением данной проблемы.
У меня работает до 2.8Msps примерно без потерь сэмплов.