[fontp 0]

В биологии сейчас это очень актуальная и востребованная задача.

Есть такой способ - цепляют к белку люминесцирующий центр и смотрят по люминесценции, куда в клетке белок пристроился.

Можно разным белкам прицепить разные центры и включать люминесценцию по очереди ( меняя длину волны возбуждения _- смотреть как белки по живой клетке

перемещаются с разрешением в нанометр.

Полагаете, это - чисто реклама?

 

Я не думаю, что здесь возможно говорить о разрешении не передёргивая понятия

А сама задача, как задача, ничего сверхъестественного

Байку про "сверхразрешение" обычно заводят, чтобы преодолеть невозможное (втихаря его переопределив, это наперсточники)

 

Ну что Вы. Форма огибающей спектра будет разной. Критерий 1/T появляется в случае, когда сигнал/шум порядка едиицы.

Если известно, что это - две синусоиды ( то есть, известна форма огибающей от каждого из сигналов в отдельности ) и сигнал/шум большой,

не надо провала в половину мощьности ( это критерий разрешения из оптики ) можно учесть существенно меньшее отличия

суммарной огибающей от огибающих каждого из сигналов в отдельности.

 

сигнал/шум не один, а должен быть тысячами и расти экспоненциально по сверхразрешению k от 1/kT.

Порядок величины 1/T не превзойти в этой жизни

 

Другое дело - позиционирование одиночного колокола. Достижима дисперсия ошибки CRB (Cramer-Rao bound)

var(f) = 12/(T*T*SNR*N*(N*N-1))

зависит линейно по SNR

[Ruslan1 17]

Доброго времени суток и всех с прошедшими праздниками! :)

 

Я никуда не пропал, тема действительна для меня актуальна но так как писать нечего то не хотел просто поднимать тему раньше чем будет что сказать. Сейчас знаю что в конце февраля у меня запланировано это бодание с новым прибором в который попробуем эту функцию влепить, тогда думаю будут и новые вопросы и новая информация. До тех пор я в-общем-то загружен другим и вряд ли вставлю что-нибудь больше чем "ага!" или "ну кто бы мог подумать!" в дискуссию. :)

 

Другое дело - позиционирование одиночного колокола. Достижима дисперсия ошибки CRB (Cramer-Rao bound)

var(f) = 12/(T*T*SNR*N*(N*N-1))

зависит линейно по SNR

 

С точки зрения моей конкретной задачи мне хочется знать, что частота синусоиды изменилась при ее изменении на величину 0.001 Hz, то есть вычисленная частота для этой новой синусоиды будет отличатся от частоты для старой синусоиды. Ничтоже сумняшеся я продолжаю называть эту определяемую разницу частот "Resolution", и переводите как сочтете нужным :)

 

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

 

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

[blackfin 32]

С точки зрения моей конкретной задачи мне хочется знать, что частота синусоиды изменилась при ее изменении на величину 0.001 Hz, то есть вычисленная частота для этой новой синусоиды будет отличатся от частоты для старой синусоиды.

А где Вы собираетесь брать опорный генератор с нестабильностью 0.001 Hz / 2*6500 Hz = 7.7e-8 ≈ 0.08 ppm?

[vallav 0]

А где Вы собираетесь брать опорный генератор с нестабильностью 0.001 Hz / 2*6500 Hz = 7.7e-8 ≈ 0.08 ppm?

 

Так ему разрешение нужно а не точность.

Подойдет генератор, который между двумя измерениями уходит не больше, чем Вами озвучено.

Если между измерениями не блее часа - такой генератор не проблема...

 

[fontp 0]

С точки зрения моей конкретной задачи мне хочется знать, что частота синусоиды изменилась при ее изменении на величину 0.001 Hz, то есть вычисленная частота для этой новой синусоиды будет отличатся от частоты для старой синусоиды. Ничтоже сумняшеся я продолжаю называть эту определяемую разницу частот "Resolution", и переводите как сочтете нужным :)

 

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

 

 

В принципе - 0.001 возможно при хорошем отношении сигнал/шум и большом числе точек. Оценка максимального правдоподобия (CRLB) - достижима, хотя её нельзя превзойти, и поэтому она всегда используется для сравнений в качестве стандарта.

 

Но играет роль не только сигнал/шум, но и помехи. Поскольку такая высокоточная оценка держится на априорном предположении, что синусоида - одна (или может быть изолирована от помех спектральными окнами с быстро убывающими хвостами)

 

Вот здесь приводятся худшие оценки CRLB для двух синусоид, находящихся рядом

http://home.uludag.edu.tr/~dilaver/web_pub..._journal_10.pdf

хотя это должно быть и интуитивно понятно, что синусовая помеха в релеевском диапазоне (отличие частот ~1/T) испортит хвостом любой алгоритм, независимо от сигнал/шум и количества измерений.

Предельная точность измерения частоты синусоиды ограничивается шумом, синусовыми же помехами и стабильностью самой обмеряемой синусоиды. Не говоря уже о стабильности сэмплинга

А где Вы собираетесь брать опорный генератор с нестабильностью 0.001 Hz / 2*6500 Hz = 7.7e-8 ≈ 0.08 ppm?

[vvs157 0]

Прочитал, но не понял. То есть у меня есть линейка с рисками 0.001 Гц, а я с ее помощью измеряю с погрешностью 0.84 Гц, так, что ли?

Как продолжение аналогии с метровой линейкой и миллиметровыми делениями. У линейки есть коэффициент теплового расширения (который на ней обычно не указывается), и если вы этой линейкой что-то измерили с точностью до мм на морозе -40 и в бане при +60, то реальная погрешность будет куда больше 1 мм.

[=GM= 0]

1.3 мм для стальной линейки. И надо учитывать коффициент линейного расширения материала объекта измерения. Если измеряете длину такой же линейки, то вгрубе разница будет пара микрон.

[Ruslan1 17]

И снова зравствуйте! :)

Всех поздравляю с прошедшими праздниками, начиная с нового нового и старого нового годов и до 8-го марта включительно :)

 

Вот и добрался я до поля мной непаханного, а именно до вычисления частоты не подсчетом переходов через нуль, а методами ЦОС.

Напомню чего лично мне хочется:

Вид сигнала- единственная затухающая синусоида (колебание струны)

Диапазон частот входного сигнала- 400-6500 Hz (ну может от 100 Hz, но сейчас не нужно).

Длительность выборки- 300 ms

Цель номер раз: определить частоту сигнала с разрешающей способностью 0.001 Hz. Все остальные параметры сигнала (амплитуда например) сейчас не интересуют.

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

 

1. Поставил (наконец-то) Матлаб,

2. оцифровал реальный сигнал через вход Sound Card. Использовал разную частоту семплирования, но всегда 4096 сэмплов. Все точки сэмплирования являются значащими (то есть синусоида присутствует и нет дополнительных шумов связанных с ее возбуждением). Входной сигнал имеет частоту в окресности 820.2 Hz

3 вычислил FFT. Ниже приведена картинка вычисленного. Четко виден основной сигнал на 820 Hz и его третья гармоника на 2460 Hz. Также видна палка на пятой гармонике (4100 Hz). палка около 4600 Hz- это шумы преобразователя питания реальной схемы, будем бороться. Что-то из шума наверняка является приветом от звуковой карты, которая и оцифровывает.

 

 

Некоторые вопросы:

1. Я специально никаких окон не задавал, Как я прочитал, матлаб по умолчанию не использует оконных функций, "размывающих" спектр основной палки. Как я понял, окно нужно выбирать в соответствии с тем, какой метод будет использоваться для уточнения частоты. Так ли это?

2. Потерялся в методах. Первый этап понятен- беру FFT и нахожу область, в которой буду делать уточнение (+/- половина максимального бина). А вот дальше туплю. Единственное что понял- математику не тяну :(. То есть из теоретической статьи написать m-файл чтобы посмотреть в работе, не смогу. Так сказать, определил пределы своих способностей :( Может быть у кого-то есть матлабовские примеры на этот счет, чтобы глянуть, или ссылки на статьи для практикующих инженеров, которые плохо и давно учили ЦОС?

Сейчас вот открыл пару десятков закладок, в том числе и электроникса, обсуждалось уже тут подобное. Может умнее стану когда прочитаю........... В-общем, определяюсь с методом и нану копать интернет на предмет доступной моему пониманию практической реализации.

 

пока что самое мне понятное, это рекомендации fontp (http://electronix.ru/forum/index.php?showtopic=80460)

Существует более перспективный подход, когда определяется максимальный бин DFT, и два скалярных произведения с комплексными экспонентами (Герцелями, если так больше нравится называть) отстоящими от частоты центрального бина на четверть ширины бина. Дальше снова проводится квадратичная интерполяция спектра мощности по 3-м точкам. В англоязычной литературе такой способ называют часто для убедительности "maxlikelihood extension", поскольку вблизи максимума любую функцию можно приблизить параболой. В принципе этот метод ничем не отличается от квадратичной интерполяции с добавлением нулей, но спектр вне точек DFT оценивается только вблизи максимального бина

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

 

[Ruslan1 17]

В-общем, я получил положительный результат

 

1. Использую FFT8192 16-битную, с фиксированной точкой, 32150 SPS

2. Окно Хэмминга

3. параболическую интерполяцию в области максимального бина.

 

результат очень замечательно болтается в диапазоне +/- 0.002 Гц :) Так что у меня 0.001 Гц resolution уже есть :)

 

Вижу одну непонятку- разница по сравнению с методом подсчета периода составляет около 0.03 Гц. Сейчас не готов сказать кто прав, но пока что думаю что метод на базе FFT ошибается.

Исходя из этого посоветуйте пожалуйста, какое из направлений более предпочтительно для улучшения результата:

1. Попробовать другое окно. Пока что пробовал без окна, Хээна и Хэмминга. Без окна жуть (на полгерца отличалось от метода измерения периода, болталось на те же пару 0.001Гц но около неправильного значения). Между Хэнном и Хэммингом разницы не заметил.

2. Попробовать перейти с 16-битных данных на 32-битные (у меня 24-битный АЦП, нормализую до максимально возможного 16-битного диапазона). Проблема- ОЗУ на FFT8192 в этом случае не зватит, только FFT4096

3. Попробовать перейти от целочисленного к плавающеточечному FFT. Опять же тогда осилю только FFT4096.

4. выбрать другой метод аппроксимации.

5. Оставить все как есть, списав на корпускулярные паразитные излучения из соседнего подпространства :) Благо нужный resolition уже достигнут, а новые хотелки в следующей итерации проверить, впаяв контроллер с бОльшей RAM.

 

Понимаю, 5-й вариант всегда предпочтительной, но может можно что-то еще.....

 

[Alex11 6]

При этих точностях 16 бит категорически не хватает, так что 32 бита или плавучка, с последней проще. Нам пришлось кое-какие внутренние переменные при счете делать double. 4096 должно хватить.

[Ruslan1 17]

При этих точностях 16 бит категорически не хватает, так что 32 бита или плавучка, с последней проще. Нам пришлось кое-какие внутренние переменные при счете делать double. 4096 должно хватить.

Большое спасибо за совет!

Попробовал подсунуть моему вычислителю 16-битный массив с идеальной синусоидой с частотой 823.150 Гц - получил результат 823.095 Гц. :(

 

Бум копать, причем просто на PC в Билдере. Рано я к железу перешел, преждевременно :(

[Ruslan1 17]

При этих точностях 16 бит категорически не хватает, так что 32 бита или плавучка, с последней проще. Нам пришлось кое-какие внутренние переменные при счете делать double. 4096 должно хватить.

Неа, вроде бы не в этом дело.

Прогнал на Матлабе, используя его встроенный fft и подсовывая идеальный синус (массив double). То есть как бы все считалась с большой точностью.

Результаты таковы (подсовываю 823.15 Гц)

Моя железяка (16-битовые целочисленные исходные данные и такие же целочисленные расчеты окна и fft) : 823.0948 Гц

Матлаб, у которого и входной массив double, и все расчеты тоже в double : 823.0944 Гц.

То есть переход от интов даже к даблам изменил результат на 0.0004 Гц.

Вывод: дело не в точности а в стратегии. Что-то у меня с окном или с методом аппроксимации делать нужно.

Может быть точность тоже вылезет, но сейчас она не виновата.

(Сама аппроксимация в float32 делается, но все расчеты до нее целочисленные)

 

PS. Однако сам слегка офигел от такой точности целочисленных вычислений. Причем выбранная частота далеко от бина.

 

[ivan219 0]

Ruslan1 вы можете выложить исходники на билдере?

И все расчёты по параболической интерполяции.

Изменено 20 марта, 2011 пользователем ivan219

[ivan219 0]

Ruslan1 Я вот что заметил. Если делать комплексный БПФ из комплексного сигнала, то спектр получается правильнее, чем, если делать из вещественного сигнала.

Пример: частота сигнала попадает ровно в середину между двумя бинами БПФ.

 

F = 127.5 при N=512

Значения при комплексном сигнале:

2122,0958687503

6366,20771054087

6366,20771054087

2122,0958687503

Как видно они зеркально расходятся от центра.

Вещественный сигнала:

2122,30561555389

6366,35751840486

6366,11782864258

2122,06590774154

Симетрия нарушена.

 

А вот как выглядит спектр при F = 2.5 N = 512

Комплексный:

2122,0958687503

6366,20771054087

6366,20771054087

2122,0958687503

Вещественный:

3031,32298384062

7073,35328007469

5787,65222982185

1632,55815496941

 

Этот метод даёт ошибку именно из за того что делается БПФ вещественного сигнала. Попробуйте с комплексным.

Изменено 20 марта, 2011 пользователем ivan219

[Ruslan1 17]

Отмечаю сразу всем, извините что не сразу.

 

0. Матлаб рулит! в который раз была доказана простая истина, что все нужно проверять на 100%, а потом уже в железо пихать :)

1. Обнаружена методическая погрешность из-за использования оконной функции. Погрешность имеет синусоидальный характер с периодом пол-бина и достигает например 0.06 Гц под Хэммингом. Вылечил использованием Гауссовского окна с альфой больше чем 3.9 - погрешность имеет тот же вид, но максимальное значение меньше чем 0.001 Гц (я стараюсь достичь такой точности).

 

2. На чистой рассчитанной синусоиде и синусоиде от генератора никаких дополнительных ошибок из-за использования целочисленного БПФ выявлено не было. Но при исследовании настоящих сигналов от датчика обнаружена разница в вычислениях в даблах Матлаба и целочисленного БПФ. Перешел на 32-битный float (беру FFT4096), разница с Матлабом стала исчезающе мала. Думаю, что реальный сигнал имеет отклонения от идеальной синусоиды, и эти искажения по-разному влияли на результат при использовании вычислений разной точности (float против short)

 

3. Насчет преобразовать в комплексный сигнал и в таком виде считать. Отличная идея, сам тож подумал что некузяво сейчас действую- заполняю Im нулями, и из результата только половина массива используется. Но так четко сформулировать не смог. :). Сейчас имею материалы как это сделать, но буду делать позже, на сегодняшний день работает, ну и ладно. Надо брать пример с Исимбаевой, которая планку чуть-чуть подымает и опять за новый мировой рекорд все заслуженные почести имеет. А прыгнула бы сразу на пол-метра выше- и все, карьере конец ;) Вот и я тоже, идеи есть, но если все сразу реализую- то завтра придется новые придумывать:)

 

4. Ошибка из-за подсовывания "псевдо-комплексного" сигнала (Im=0) комплексному БПФ: странно как-то. Постараюсь проверить, но в теории какая разница, ну легли у меня векторы на ось, почему это криминал?

 

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

 

Озадачило то, что все равно вижу разницу при измерении реального датчика (не генератора синусоиды) порядка 0.02 Гц между методом подсчета пересечений нуля и на базе FFT. Но проблема в том, что между измерениями есть промежуток 2 секунды, может просто датчик меняет свои показания. Проверил на Матлабе возможную дополнительную погрешность если входной сигнал не идеальный синус, а затухающая синусоида- получил очень малые отклонения (до 0.001Гц) при скорости затухания в десятки раз большей чем на моем датчике.

Может быть кто-то может помочь с проверкой? Могу предоставить массив 24-битных семплов c АЦП, Хочу узнать точную частоту, которая получается при вычислении (любым Вашим методом) . Если будет отличаться от того что я насчитал- значит все-таки ошибка у меня в методе.