Здравствуйте!

Помогите пожалуйста, спать не могу, все думаю......

Накопал в интернете, не понимаю как сделано. а хочется понять и применить методу.

короткая листовка: http://www.campbellsci.com/documents/produ...es/b_avw200.pdf
полный мануал: http://www.campbellsci.com/documents/manuals/avw200.pdf


Прибор определяет частоту синусоидального сигнала, который может быть в диапазоне 100 - 6500 Гц, с разрешением 0.001Гц и точностью 0.013% от измеренной величины. Указано что используется 4096-точечное FFT. Интересно, что еще кроме упомянутого FFT применяется для получения такого разрешения по частоте?

Как это может быть сделано? Взяли FFT, определили частоту грубо, потом детально обнюхали область вокруг найденной частоты, получили частоту точно?


Что абсолютно точно известно (это физика процесса):
1. Полезный сигнал- честная одночастотная синусоида (дернули струну, потом ищут ее частоту собственных колебаний. то есть это затухающее колебание механической струны, которое преобразовано в электрический сигнал). Шумы возможны.
2. Общая длительность исследуемого сигнала доли секунды (ну пусть 300 ms)


Что говорит изготовитель (интересный текст и какие-то картинки начинаются со страницы 85, еще есть спецификация на странице 9):
1. The AVW200 uses an audio A/D for capturing the sensor’s signal. The number of samples acquired in this period is 4096 points. A Fast Fourier Transform (FFT) algorithm is used to create a frequency spectrum.
2. the spectral analysis gives improved frequency resolution (0.001 Hz rms) during quiet conditions.
3. the AVW200 Fourier transforms the recorded response and analyzes the resulting spectrum to determine the wire’s resonant frequency. This analysis also provides diagnostic information indicating the quality of the resonant-frequency measurement.
4. Measurement resolution: 0.001 (Hz RMS)
5. Accuracy basic: +/- 0.013% of reading.
6. DF measurement time between 1.6 to 1.8 second
7. SYSTEM: PROCESSOR: Hitachi H8S 2324 (16-bit CPU with 32 bit internal core), MEMORY: Either 128 or 512 kbytes of SRAM; 2 Mbyte of OS Flash
8. частоту проца не знаю, но пишут про дополнительный ток потребления во время измерения 27mA@12V. Судя по даташиту проца получается что используется он по максимуму, 25 МГц тактовой.


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

может быть применен какой-н сплайн к спектральным составляющим fft
под п.2 можно понимать что угодно, хоть нейронные сети :D ))

За 300 мс можно измерить "методом захвата" с точностью 0.001 Гц, поищите, была такая тема в АВР МК.. и спите спокойно.

Это переходы через ноль считать (zero-crossing counting method)? Не, неинтересно. Так работало старое поколение таких измерителей, их ахиллесова пята это шумы из-за которых лишние переходы. Кстати там разрешение на порядок меньше всегда анонсировалось (0.01Гц).

А про точность 0.001 методом захвата- это вы я так думаю погорячились. Ну во-первых нужно переименовать "точность" в "разрешение" (ну или говорить о стабильности 0.15ppm, что за АВР такой?), а во-вторых при 6.5 кГц это будет 23 пикосекунды разница в длительности периода. 44 наносекунды за 1950 периодов (300мс и 6500Гц). это 22 МГц тактирование таймера чтобы один такт набежал. Ладно, тут наверное можно и поверить, уж это AVR наверное может.



Спасибо, попробую погуглить в этом направлении. Но опять же смущает заявленная у них точность 0.013%. Я думал что любая апроксимация может увеличиить только разрешающую смособность, но не гарантировать точность. А у них получается что точность на уровне 1/7700, то есть почти в 4 раза выше чем если за точность брать расположение палок исходного FFT.
Но тут наверное хитрость в том что точно известна форма искомой функции- чистый синус и это может сыграть.
Я вот еще подумал, если после FFT я знаю узкий диапазон и делаю цифровую фильтрацию, отсекая все выше и ниже, может это помочь при дальнейшей обработке? но вот непонятно что за бработка такая....
Спасибо за подсказку, завтра про сплайны почитаю.

А насчет нейронных сетей- не та это техника. Это полевой прибор, не для лабораторий. Думаю все сделано по максимуму дубово и без новомодностей.

Тут, действительно, фишка в том, что сигнал - чистый синус. Здесь нет спектрального разрешения 0.001Гц в смысле, что две линии на таком расстоянии друг от друга отличить будет невозможно. Имеется в виду, что младшая значащая цифра прибора 0.001Гц. Мы сделали прибор для измерения синуса в сети, так там разрешение и точность еще выше. Использовано понятно что. На входной сигнал накладывается оконная функция. Что они использовали - не знаю, а у нас - Гаусс. Далее, форма линии в спектре будет совпадать с формой оконной функции. Подгоняем коэффициенты так, чтобы функция наилучшим образом легла на результат Фурье-преобразования. Из коэффициентов считаем точную величину частоты. Получается очень хорошо.

Именно так..

1. Взяли FFT4096 от 0 Hz до 16384 Hz и нашли f_min и f_max,
причем, ∆f = f_max - f_min = 16384 Hz / 4096 = 4 Hz;
2. Умножили входной сигнал f_s на комплексную экспоненту с частотой -f_min.
3. Отфильтровали сигнал в полосе 0 Hz до 4 Hz.
4. Сделали децимацию в 4096 раз.
5. Добавли четыре полученных семпла в буфер FIFO на 4096 точек. Четыре первых семпла из буфера FIFO удаляем.
6. Взяли FFT4096 от 0 Hz до 4 Hz и нашли f_s = f_min+f_bin,
причем, ∆f_s = 4 Hz / 4096 = 0.0009765625 Hz;

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

p.s. про нейронные сети шутка такая, их вставляют куда не попадя в рекламу, как бананотехнологии

если вы хотите получить точность 0.001 Гц при использовании fft при неизвестном сигнале, то вы должны наблюдать сигнал 1000 сек. Поэтому по-любому используют дополнительную информацию для повышения точности. Я так понимаю, что исследуется колебание струны , а значит форма сигнала известна (затухающая синусоида), неизвестны лишь несколько параметров, частота, нач фаза, коэффициент затухания или что-то в этом роде. Остается произвести анализ параметров вашего сигнала. Для этого есть куча вариантов, например методом наименьших квадратов путем оптимизации подобрать эти параметры. Если нужна только частота, то можно амплитуду ограничить и ставить PLL, нейронные сети предлагали, в общем поле для фантазии.


Я и говорю надо ждать 1000 сек не думаю что это возможно. Ибо сигнал затухает и нет голой синусоиды а есть некий горбик.

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

Благодарствую!
у меня путем несложных калькуляторных вычислений на уровне интуиции изначально крутилось что-то валенковское типа "два раза применили FFT4096 к входному реалу, уж больно в цифры хорошо попадает". Но интуиция на пустом месте дело бесполезное.... буду копать в указанном направлении. Хотя какое тут направление, готовый алгоритм, сразу кодировать можно


2 All:

Благодарю всех за дельные советы! Вы мне показали что такое возможно и дали даже больше одного варианта как это сделать.

В-общем сказал начальству что умные люди сказали что оно не бред и реально, по крайней мере есть методы
Для проверки даже заложим в прототип (в котором по-старинке считать собирался) какую-никакую вычислительную мощность в виде dsPIC и на всякий случай внешнего RAM (хотя надеюсь внутренних 16К хватит) и вперед. Ну и дополнительные три недели только на проверку и реализацию этого алгоритм затребовал

Разумеется сначала попробую на компьютере прогнать снятый с АЦП сигнал (такие частоты и саундбластер возьмет), потом уже в железо перекину.

Еще раз спасибо!

Уверен что появятся новые вопросы- так я ж молчать не буду.... Но это уже после пьянки празников. У нас тут обычно и католическое Рождество мимо не проходит, чувствуется близость запада.....

Нет, не переходы. Считается целое число периодов входной частоты. Заполнение - 20 МГц, отсюда легко вычисляется погрешность.


1) Не пойму, причём здесь стабильность 0.15ppm и АВР? Подайте на АВР стабильный клок, всех делов. Если есть нестабильность опоры, значит, как по теории, складывайте погрешность метода и погрешность опоры, другого не дано.

2) Отчасти вы правы, 0.001 Гц для метода захвата в вашем случае это не точность и не разрешение, это погрешность измерения (максимальная ошибка), т.е. ошибка результата измерения будет не хуже. (Не хотел говорить, а то вы спать не будете, но всё-таки скажу: на таком интервале я могу легко получить погрешность 0.0001 Гц).

Не, теперь уже спать буду. много вариантов уже предложено, ну хоть один-то сработает
Если нет- то конечно вернемся в прошлый век, периоды считать.
Я сам против использования тяжелой артиллерии при отстреливании мелких пернатых, но тут мне кажется ЦОС на своем месте будет, да и некоторые дополнительные вкусности дает в виде измерения амплтитуды и отношения сигнал/шум. Это может предсказать отказ до выхода из строя (датчики часто вмурованы в стены во время строительства или болтаются на дне морском или еще где в глубоких скважинах и если их замена может прогнозироваться и осуществляться планово а не аварийно- это прямо праздник!). Лично я хочу применить ЦОС, но не знаю как (но теперь уже имею конкретные пути для раскопок .


О! точно все-таки юзают они FFT! они там еще в этом приборе и амплитудное значение и SNR считают! и все это на базе одной выборки за время одного измерения (каждые 2 секунды новое значение считают, а вроде бы чаще возбуждать струну не рекомендуется).

два раза fft нельзя для коротких сигналов, fft - дает усредненные значение на интервале. Получите усреднение на 1024 секундах, че там вылезет не угадать
fft + wavelet,
fft + pll
или валенком - fft + уточненный максимум по соседним спектр.компонентам

Точность измерения частота одиночной синусоиды за любое время больше периода ограничена только отношением сигнал/шум (согласно оценки максимального правдоподобия Крамера_Рао) и собственно стабильностью самой частоты. Предельная точность реально достигается с помощью интерполяции спектра. Такшта нивапрос))
ТАК!

Здесь исчерпывающий отчет по теме музицирующих программистов из Стэнфордского университета. В форуме обсуждалось неоднократно
https://ccrma.stanford.edu/STANM/stanm/node3.html
Stanford University Department of Music (STAN-M)
Technical Reports On-Line:
Abe, M., and J. O. Smith. 2004c.
Design Criteria for the Quadratically Interpolated FFT Method (I): Bias due to Interpolation.
Tech. rept. STAN-M-114. Stanford University, Department of Music.
Abe, M., and J. O. Smith. 2004d.
Design Criteria for the Quadratically Interpolated FFT Method (II): Bias due to Interfering Components.
Tech. rept. STAN-M-115. Stanford University, Department of Music.
Abe, M., and J. O. Smith. 2004e.
Design Criteria for the Quadratically Interpolated FFT Method (III): Bias due to Amplitude and Frequency Modulation.
Tech. rept. STAN-M-116. Stanford University, Department of Music.
Abe, M., and J. O. Smith. 2004a.
AM/FM Rate Estimation and Bias Correction for Time-Varying Sinusoidal Modeling.
Tech. rept. STAN-M-118. Stanford University, Department of Music.
Abe, M., and J. O. Smith. 2004b.
CQIFFT: Correcting Bias in a Sinusoidal Parameter Estimator based on Quadratic Interpolation of FFT Magnitude Peaks.
Tech. rept. STAN-M-117. Stanford University, Department of Music.

нету у меня много секунд, у меня сигнал затухающий. Я так понял что blackfin имел в виду обработку однажды ссобранных данных, два fft но второе применяется к той же изначальной выборке но уже обработанной. То есть на все-про все я имею одну исходную выборку из которой все и решаю.

Да, наверное для красоты и ясности нужно было картинку сигнала прямо тут привести. Это после того как струну дернули она так вибрирует. Сколько миллисекунд используется для выборки не знаю но получается что-то никак не дольше 300ms, даже может 200. Это хорошо ложится на 4096 точек и обычный средний 24-битный АЦП
Изготовители написали, что это сигнал снятый с "Good Sensor with a Wider Range (200 to 6500 Hz)". Это из того документа, на который я ссылался в самом первом письме.

Нажмите для просмотра прикрепленного файла



Спасибо! не пройду мимо! На подобное и надеюсь, что вопрос не нов и где-то уже все есть...... Но уж никак не думал, что "Department of Music" интересуется подобными точностями, даже если в универе. наверное математикам в общаге спать мешали

Не новый. Всё придумано до нас и даже до этих американцев и прочих немцев из соседней темы.
Я, например, использовал квадратичную интерполяцию (для определения частоты несущей модема) за 5 лет, до этой публикации музыкальных американцев.
И тогда это уже было не ново. Теория была разработана в 80-х годах. Rife&Burstin были первопроходимцами в этом вопросе ещё в 1974-м году
Пионерская статья Rife&Burstin
Но в отчетах по ссылке очень много полезных практических экспериментальных деталей. Судя по отчетам там на факультете хорошие музыкальные теоретики теормуза

Ruslan1, не хочется встревать в перепалку, но вот в посте #10 было два вопроса в ответ на ваш пост #4, вы их как бы не заметили и обошли, но зачем-то ответили только на фразу в скобках, которая совсем не главная.

И ещё вот, интересно узнать про разрешение и точность для обычной металлической линейки длиной 1м. Какое у неё разрешение и точность?

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

0. "Нет, не переходы. Считается целое число периодов входной частоты". считать не переходы а целые периоды. Если имеется в виду вычисление методом подсчета однотипных (только задних или только передних) фронтов откомпарированного сигнала- Да, согласен, в комментариях не нуждается. Это точнее засчет исключения погрешность от сдвига на постоянную составляющую. Но не более. Сейчас попробую накопать тему на которую вы ссылались, может имелось в виду что-то похитрее. Каюсь, не искал и соответственно не смотрел. Думал что тривиальный счетчик, но сейчас посмотрю чтобы не оставить без внимания может быть хороший метод.

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

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

Про линейку не забудьте. Заодно ответьте себе, можно ли с помощью такой миллиметровой линейки измерить длину в 0.08 мм? Конкретно я знаю, что можно.

Вы не поняли. Человек как бы намекает, что то что Вас интересует, в названии темы ошибочно называется разрешением.
Разрешение - фундаментальное понятие в физике, относится оно к различению двух соседних объектов.
А при определении местоположения одного объекта уместно говорить о точности или погрешности, чтобы не вызывать путаницу,
называя одним словом разные вещи, а разными словами одно и то же) Да ещё и не теми как принято по конвенции


Чтобы обеспечить РАЗРЕШЕНИЕ в 0.001 гц нужно как минимум 1000 сек, как было уже сказано выше. Разрешение ограничено принципом неопределённости,
в то время как точность измерения (дисперсия ошибки) только границей Крамера-Рао. Причем для несмещенной оценки систематическая ошибка отсутствует и это возможно.

Теперь понял, спасибо. А вот намеков не понял, отупел совсем. Это был мой вольный перевод иностранного слова "resolution". А под точностью я раньше подразумевал "accuracy" и не удивлялся тому что это разные вещи, причем между ними всегда есть определенное соотношение.

В-общем, просто сообщаю о задаче что знаю, не впадая в переводческий грех:

Measurement Resolution : 0.001 (Hz RMS)
Accuracy Basic: ±0.013% of reading (осмелюсь перевести это в герцы, зная что максимальная читаемая величина 6500Гц : ±0.84 Гц )

Measurement Resolution - это разрешение шкалы (линейки, а не прибора)= абсолютная точность.
Accuracy Basic- это относительная точность или погрешность.

Я уже признал в соседнем посте что с улавливанием намеков у меня туго. Понял. Можно. Вы знаете как это можно. Я ж разве против?


Прочитал, но не понял. То есть у меня есть линейка с рисками 0.001 Гц, а я с ее помощью измеряю с погрешностью 0.84 Гц, так, что ли?
Извиняюсь если резко звучит, но мне это непонятно. А хочется понять, чтобы не ставить перед собой недостижимых задач, которые в дейчствительности никто не решал.

Иду читать про погрешности и линейки.

PS вот так и рушаться стереотипы..... жил себе не тужил.....


Так это они что, имеют в виду что на полученной ими спектрограмме палка установлена с точностью до 0.001 Гц (расстояние между соседними палками не сообщается, просто координаты этой палки известны с указанной точностью), а вот расстояние между показываемым значением и реальной частотой - не более чем 0.84 Гц ?
Это правильно? Или опять мусор в голове?

Типа того, сказано же ±0.013% от измеряемой величины. От 65000 - 0.84, но от 6500 - 0.084


Нет не правильно.
Поскольку первая величина RMS - это и есть случайная ошибка метода измерения. Они риски так поставили, как даёт метод (при некоторых условиях, поскольку
существует фундаментальный принцип Крамера-Рао и этот предел нельзя превзойти, а он зависит от отношения сигнал/шум, который не всегда возможно контролировать)
Вторая величина - систематическая величина смещения. То что нет фундаментальных принципов, запрещающих делать несмещенную оценку, не означает,
что оценка конкретного метода несмещенная - она смещенная из-за ограничений любого метода. Из-за принятых упрощений, или недостаточной информации или
вычислительно сложно. Делайте всегда оценку максимального правдоподобия честно и всё будет несмещённо, так не хотят нелинейщины ))

Например, метод квадратичной интерполяции, про который я говорил требует применения окон и добавления очень большого числа нулей к ДПФ, чтобы убрать смещение
в ассимптотике )

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

fontp, ну что же вы? Лишили парня сна на целую неделю :-)

Если основная частота одна (либо несколько, но разнесены достаточно широко), можно сделать проще, без второго FFT.

Мы делали прибор для метрологических измерений на сетевой частоте (50 +/- 2Гц). Выборка - 300 мсек - примерно 15 периодов, , самплинг 25600 Гц, 24-битное АЦП, для обрубания концов - использовали гауссовское окно. В результате по одной выборке частота определялась до, примерно, 10^-4 - 10^-5 Гц, это реальные полученные цифры.

Гауссовское окно для моночастоты дает опять же гауссовское частотное распределение с затуханием порядка 80 дБ. При данной выборке ширина порядка 30-40 Гц (точно уже просто не помню). Этот гаусс по МНК аппроксимировался гауссом для нахождения максимума по ближайшим 6-ти частотным точкам (влево и вправо). В результате получали приведенную выше точность даже при высоком уровне помех.

Для расширения частотного диапазона 36 - 12000 Гц делался дополнительный шаг - сначала находилась частотная компонента с максимальной амплитудой, затем по тому же алгоритму считали в ее окрестности. Результат аналогичен.

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

Вы что то путаете. Принцип неопределенности не имеет касательства к координатам двух разных объектов.
Он связывает или координаты и импульс или энергию и время для данного объекта.

Если известно, что в сигнале две синусоиды, то разрещение ограничено величиной сигнал/шум и может быть много меньше
1/T.

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




Разрешение - это возможность отличить один объект от другого. Относительная величина,
сравнение друг с другом.
В данном случае - если частота синусоиды изменилась - какое изменение будет замечено.
Точность - это возможность определить, чему именно величина равна - абсолютная величина.
сравнение с эталоном.
Обычно разрешение бывает намного больше точности.

1. Неоднозначность локализации любой функции связана с неоднозначностью локализации её фурье-образа принципом неопределенности-
dT * dF > 1
Откуда и берется в физике то о чём Вы говорите про импульс)) Принцип неопределённости появился не в квантовой механике, а уже в Фурье-оптике. Ну а теперь прикиньте как Вы будете "разрешать координаты" (частоты в данной задаче) двух объектов, если они в принципе размазаны на 1/T. Другое дело - измерить, если объект был один. Да хотя бы взять центр масс кривой его рассеяния (кружка рассеянияния в оптике-так и делают в звездных координатных датчиках). А ещё лучше подогнать под известный образ кривой с помощью интерполяции и достигнуть предела Крамера-Рао.



2. Нет. Если известно, что в сигнале одновременно присутствуют две синусоиды - Вы не сможете их отличить если они рядом, ближе чем 1/T и этот факт от отношения сигнал/шум практически не зависит. Вы можете этого не знать, но доказано, что если в сигнале 2 (несколько) синусоид и они расположены достаточно далеко друг от друга, то предельная точность измерения действительно зависит от отношения сигнал/шум в соответствии с той же предельной оценкой Крамера-Рао для одиночной синусоиды. Но когда они сближаются поближе к 1/Т возможность одновременной оценки их частоты резко теряется. Существует формула Крамера-Рао для предельной оценки по максимуму правдоподобия для ДВУХ КОМПЛЕКСНЫХ СИНУСОИД. Об этой формуле редко вспоминают. Просто она практически бесполезна ввиду своей громоздкости. Но теоретическую ценность она представляет, поскольку показывает как летят к чертовой матери любые попытки оценки частот двух синусоид при их сближении по частоте... Я дам Вам ссылку - если Вам, конечно, захочется))



3.Вы можете измерить положение одиночной сингулярности с очень высокой точностью, ограниченой только критерием Крамера-Рао оценки максимального правдоподобия. Не правильно называть это сверхразрешением, поскольку используется тот априорный факт, что объект был один,
а значит ничего "разрешать" не пришлось. Сверхразрешением обычно называют несколько другую вещь -повышение разрешения через решение некорретктной обратной задачи. Причем обычно она не особо решается при сколько нибудь разумных SNR.

Но Вы не можете отличить 2 близкорасположеных объекта предъявляемых одновременно. (В этом как раз и есть смысл "отличения светящихся по заказу") По закaзу оно то конечно, поскольку это одиночные синусоиды (или звезды), предъявляемые не одновременно, а последовательно. К разрешению в смысле оптики это не имеет никакого отношения. Такие произвольные толкования годятся только для рекламных буклетов. Разрешить - это именно отличить на изображении два объекта, предъявляемых одновременно.
Это значительно сложнее чем измерить положение одного объекта, поскольку их изображения размазаны и накладываются.

Ну что Вы. Форма огибающей спектра будет разной. Критерий 1/T появляется в случае, когда сигнал/шум порядка едиицы.
Если известно, что это - две синусоиды ( то есть, известна форма огибающей от каждого из сигналов в отдельности ) и сигнал/шум большой,
не надо провала в половину мощьности ( это критерий разрешения из оптики ) можно учесть существенно меньшее отличия
суммарной огибающей от огибающих каждого из сигналов в отдельности.





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

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



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

Другое дело - позиционирование одиночного колокола. Достижима дисперсия ошибки CRB (Cramer-Rao bound)
var(f) = 12/(T*T*SNR*N*(N*N-1))
зависит линейно по SNR

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

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



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

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

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

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

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

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

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

Вот здесь приводятся худшие оценки CRLB для двух синусоид, находящихся рядом
http://home.uludag.edu.tr/~dilaver/web_pub..._journal_10.pdf
хотя это должно быть и интуитивно понятно, что синусовая помеха в релеевском диапазоне (отличие частот ~1/T) испортит хвостом любой алгоритм, независимо от сигнал/шум и количества измерений.
Предельная точность измерения частоты синусоиды ограничивается шумом, синусовыми же помехами и стабильностью самой обмеряемой синусоиды. Не говоря уже о стабильности сэмплинга

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

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

И снова зравствуйте!
Всех поздравляю с прошедшими праздниками, начиная с нового нового и старого нового годов и до 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)

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

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-й вариант всегда предпочтительной, но может можно что-то еще.....

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

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

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

Неа, вроде бы не в этом дело.
Прогнал на Матлабе, используя его встроенный fft и подсовывая идеальный синус (массив double). То есть как бы все считалась с большой точностью.
Результаты таковы (подсовываю 823.15 Гц)
Моя железяка (16-битовые целочисленные исходные данные и такие же целочисленные расчеты окна и fft) : 823.0948 Гц
Матлаб, у которого и входной массив double, и все расчеты тоже в double : 823.0944 Гц.
То есть переход от интов даже к даблам изменил результат на 0.0004 Гц.
Вывод: дело не в точности а в стратегии. Что-то у меня с окном или с методом аппроксимации делать нужно.
Может быть точность тоже вылезет, но сейчас она не виновата.
(Сама аппроксимация в float32 делается, но все расчеты до нее целочисленные)

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

Ruslan1 вы можете выложить исходники на билдере?
И все расчёты по параболической интерполяции.

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

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

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

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 АЦП, Хочу узнать точную частоту, которая получается при вычислении (любым Вашим методом) . Если будет отличаться от того что я насчитал- значит все-таки ошибка у меня в методе.

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

Ага, согласен. и у меня то же самое, даже график приводить не буду- один в один. конечно не синусоида, под "синусоидальностью" я подразумевал форму, издалека похожую на синус Ну и период равен бину, а не пол-бина, как в моем сообщении, это я куда-то не туда посмотрел.

Но вопрос остается: можно ли изменением обработки приведенного массива изменить величину вычисляемой частоты?
Если кто-то посчитает по своей методе и скажет мне результат, буду очень сильно благодарен!

Вот привожу файлик исходных семплов АЦП (странно, файл с расширением csv загружать запрещено, пришлось в txt переобозвать).
Нажмите для просмотра прикрепленного файла

Файл имеет простой текстовый csv-формат, который хоть в ексель хоть в матлаб подсовывается с пол-пинка. на всякий случай привожу код, как его в матлаб прочитать:


Мои результаты вычислений:
1. максимальный бин (индексы начинаются с единицы, как принято в Матлабе):
n_bin_max = 109

2. расстояние от максимального бина:
d = -0.3976

3. собственно координаты в бинах :
bin = n_bin_max + d = 108.6024

4. Частота при указанных выше 31250 SPS и 4096 точках:
freq2 = (bin-1) * Fs / (FFT_N_SAMPLES) = 820.9413

Согласно матлабу, результат может изменится на величину до -0.02 Гц при изменении коэффициента альфа в Гауссовском окне (от 4 до 10). Но в любом случае не достигается величина, показываемая методом подсчета пересечений через ноль (разница до -0.09 Гц от FFT-based метода)

Кстати, у меня так и получается, что "результат всегда находится ближе к центральному бину чем реальное значение". Насчет "всегда" не знаю, но вот та же петрушка.....

Поанализировал я тут Ваш файлик. Очень тяжело искать черную кошку в темной комнате, особенно, еслиее там нет. Так вот. Мало того, что это не синусоида, мало того, что зашумлено, но частота сигнала изменяется на протяжении этого семпла. Так о чем может идти речь? На 1/4 семпла от начала частота 820.9137, на 1/2 - 820.9203, на 3/4 - 820.9486. За последний знак не поручусь, т.к. частота плывет непрерывно и это тоже усредненные данные. Интегрально мой алгоритм дает частоту 820.9370. И это естественно, т.к. окно отрезает начало и конец, затем Фурье дает некоторое интегральное значение.
А откуда значение 820.9413, которое хотелось получить?

Огромное спасибо! Вы меня обнадежили! Значит, я где-то напортачил с методом, если из тех же данных получаю не то. Как я говорил, альтернативный метод подсчета (измерение длительности периода с усреднением)дает меньшую величину, чем у меня. Полученная Вами цифра мне очень нравится , пожалуйста, помогите мне ее достигнуть!

Насчет несинусоидальности. Для меня струнные датчики новая область, но везде пишут что это именно синус на выходе. И датчик у меня честный промышленный, не самоделка. значит, что-то с условиями я не соблюдаю у себя на столе. Ну или как вариант- моделируется это колебание чем-то(дом шатается, я уже такое проходил с сейсмосенсорами, 7-й этаж уже болтается вполне предсказуемо), а датчик чувствительный. Это же относительно изменения частоты- может это так и должно быть, реально меняется показание.
Искажение в усилительных каскадах прибора- может быть, конечно, но вот синус с генератора нормально через них проходит....
Но в данный момент очень хочется попробовать доползти до Вашей цифры.

Как я уже писал, для меня ЦОС это новая область. Я уже не говорю о том, что Матлаб впервые установил дней десять назад.
До всего дохожу в рабочем порядке, читая форумы и интернет-ресурсы. В результате родился этот файл Нажмите для просмотра прикрепленного файла (нужно переименовать в *.m). ниже привожу его же.

он и возвращает указанное значение: Frequency_Response = 820.9413



Самым слабым местом мне кажется оконная функция, хотя может еще где напортачил.


Сейчас нашел файлик, который я оцифровывал через саундкарту с этого же датчика в тех же условиях, но с другим аналоговым трактом и другим АЦП. Все-таки больше похоже на синус. Нажмите для просмотра прикрепленного файла
Файл читается в Матлаб после распаковки командой
Значит я не только в математике напортачил..... Это я конечно завтра начну рыть....
Но вопрос остается открытым- как из того недосинуса получить частоту меньше чем у меня на полной выборке, неправильно выбран метод или более конкретные ошибки при реализации.....

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