как это сделано: FFT4096, 100-6500Hz, freq.resolution 0.001Hz, как достигнуто такое разрешение по частоте? Опции

[Ruslan1]

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

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

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

короткая листовка: 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 лет назад в институте так с тех пор ни разу ничего толком кроме чужой математики для ЦОС не применял. В-общем валенок и может быть это нужно в раздел для начинающих перетащить если что-то понятное любому ЦОС-нику спрашиваю.

[Fast]

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

[=GM=]

.. спать не могу, все думаю..

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

[Ruslan1]

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

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

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

может быть применен какой-н сплайн к спектральным составляющим fft

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

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

[Alex11]

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

[blackfin]

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

Именно так..

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;

Сообщение отредактировал blackfin - Dec 22 2010, 03:39

[Fast]

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

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

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

[bahurin]

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

3. Взяли FFT4096 от 0 Hz до 4 Hz и нашли fs = fmin+fbin,
причем, ∆fs = 4 Hz / 4096 = 0.0009765625 Hz;

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

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

Сообщение отредактировал bahurin - Dec 22 2010, 02:45

[Ruslan1]

Именно так..

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;

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


2 All:

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

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

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

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

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

[=GM=]

Это переходы через ноль считать (zero-crossing counting method)?

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

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

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

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

[Ruslan1]

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

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


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

[Fast]

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

[fontp]

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

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

ТАК!

Здесь исчерпывающий отчет по теме музицирующих программистов из Стэнфордского университета. В форуме обсуждалось неоднократно
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.

[Ruslan1]

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

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

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

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

Здесь исчерпывающий отчет по теме музицирующих программистов из Стэнфордского университета. В форуме обсуждалось неоднократно
https://ccrma.stanford.edu/STANM/stanm/node3.html

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

[fontp]

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

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

[=GM=]

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

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

[Ruslan1]

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

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

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

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

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

[=GM=]

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

[fontp]

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

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


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

[Ruslan1]

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


Чтобы обеспечить РАЗРЕШЕНИЕ в 0.001 гц нужно как минимум 1000 сек, как было уже сказано выше

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

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

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

[fontp]

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

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

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

[Ruslan1]

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

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

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

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

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

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

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

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

[fontp]

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

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

Это правильно?

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

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

[Ruslan1]

Нет не правильно.

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

[=GM=]

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

[rudy_b]

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

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

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

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

[Alexey Lukin]

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

[vallav]

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


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

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

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

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

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

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

[fontp]

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

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

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

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

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

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

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

[vallav]

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

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

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

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

[fontp]

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

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

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

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

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

[Ruslan1]

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

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

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

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

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

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

[blackfin]

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

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

[vallav]

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

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

[fontp]

С точки зрения моей конкретной задачи мне хочется знать, что частота синусоиды изменилась при ее изменении на величину 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.001 Гц, а я с ее помощью измеряю с погрешностью 0.84 Гц, так, что ли?

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

[=GM=]

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

[Ruslan1]

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

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

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]

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

[Ruslan1]

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

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

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

[Ruslan1]

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

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

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

[ivan219]

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

Сообщение отредактировал ivan219 - Mar 20 2011, 17:49

[ivan219]

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

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

Сообщение отредактировал ivan219 - Mar 20 2011, 19:33

[Ruslan1]

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

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

[ivan219]

1. Обнаружена методическая погрешность из-за использования оконной функции. Погрешность имеет синусоидальный характер с периодом пол-бина

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

Сообщение отредактировал ivan219 - Mar 25 2011, 09:01

Эскизы прикрепленных изображений

 

[Ruslan1]

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

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

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

Вот привожу файлик исходных семплов АЦП (странно, файл с расширением csv загружать запрещено, пришлось в txt переобозвать).
 ADC.TXT ( 33.43 килобайт ) Кол-во скачиваний: 119


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

Код

Fs = 31250
FFT_N_SAMPLES = 4096
%% input data
% read data from file
filename = 'ADC.TXT';
fid=fopen(filename','rt');
data_in = fscanf(fid, '%d,');
fclose(fid);

Мои результаты вычислений:
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 метода)

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

[Alex11]

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

[Ruslan1]

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

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

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

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

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

Код

FFT_N_SAMPLES = 4096;
Fs =31250;
alpha = 3.9
    
%% input data
% read data from file
filename = 'ADC.CSV';
fid=fopen(filename','rt');
data_in = fscanf(fid, '%d,');
fclose(fid);
%% windowing
for i=1:FFT_N_SAMPLES
        %Gaussian:
        tmp_n = abs(i - (FFT_N_SAMPLES)/2);
        tmp = alpha * tmp_n / (FFT_N_SAMPLES/2);
        coeff = exp(-0.5 * tmp*tmp);
        data_win(i) = data_in(i) * coeff;
end
%% fft
data_fft = fft (data_win,FFT_N_SAMPLES);
%% magnitude
for i=1:FFT_N_SAMPLES/2
    data_fft_abs(i) = 2.0 * abs(data_fft(i));
end
%% maximum searching
max = data_fft_abs(1);
n_bin_max = 1;
for i=1:FFT_N_SAMPLES/2
    if (max < data_fft_abs(i))
        max = data_fft_abs(i);
        n_bin_max = i;
    end
end
%% Gaussian interpolation
            m1 = data_fft_abs(n_bin_max-1);
            m2 = data_fft_abs(n_bin_max);
            m3 = data_fft_abs(n_bin_max+1);
            top_part = log (m3/m1);
            bot_part = 2*log(m2*m2/(m3*m1));
            d = top_part/bot_part;
% f = ( i + d ) * ( sr / n )
            n_bin_response = n_bin_max + d;
            Frequency_Response = (n_bin_response-1) * Fs / (FFT_N_SAMPLES)

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

Поанализировал я тут Ваш файлик. .... Мало того, что это не синусоида, мало того, что зашумлено..........

Сейчас нашел файлик, который я оцифровывал через саундкарту с этого же датчика в тех же условиях, но с другим аналоговым трактом и другим АЦП. Все-таки больше похоже на синус.  in_data32000SPS_4096pnt.zip ( 8.08 килобайт ) Кол-во скачиваний: 73

Файл читается в Матлаб после распаковки командой

[data,time] = daqread('in_data32000SPS_4096pnt.daq');

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

[Alex11]

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

[ivan219]

Поможет БИХ фильтр. Если есть вычислительные ресурсы, то хороший фильтр избавит от гармоник. Единственно, что плохо так это вычисление в несколько этапов.

[Alex11]

Фильтр не поможет, т.к. проблема не в гармониках - они и так хорошо разделяются Фурье-преобразованием, а в том, что частота плывет. И в зависимости от того, в какую часть семпла мы смотрим лучше (где наложено окно), тот ответ мы и получаем.
Теперь про результаты анализа файла с 32к оцифровки. Тяжело плавать в серной кислоте с оторванными ногами. У меня матлаб не стоит, так что пришлось потрошить файл вручную. Где-то я ошибся на один семпл и у меня прошел сбой по частоте. Но в целом этот файл гораздо лучше последнего, частота болтается не больше, чем 2 единицы второго знака после запятой, хотя амплитуда маленькая и шумов довольно много. Если сохраните файл в текстовом виде, могу посчитать еще раз поточнее. Первая прикидка частоты - 822.557.

[Ruslan1]

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

Да конечно пожалуйста. Сконвертировал из *.daq и добавил в архив и текстовый файл, расширение "csv".  in_data32000SPS_4096pnt.zip ( 22.54 килобайт ) Кол-во скачиваний: 48


И опять засада: у меня 822.5698 Hz

четверти от выборки (по 1024 точек):
1: 822.5570
2: 822.5868
3: 822.5385
4: 822.5531
половинки от выборки (по 2048 точек):
1: 822.5842
2: 822.5649
Полная выборка (4096 точек): 822.5698

Получается, что измеряемая частота сильно плывет? такое теоретически возможно. Но как-то не вяжется с тем что я вижу глазами: я последовательно применяю два метода (подсчет периодов методом пересечение нуля и FFT-based: первый-второй-первый-второй...), и вычисленная каждым из них величина колеблется около своего значения частоты, причем FFT-baised больше. И разность того же порядка как разница между Вашей и моей цифрами.

[Alex11]

Вот уточненные результаты: средняя частота 822.5666. Если смотреть по частям (я, правда, делал несколько иначе - оставлял FFT на 4096, но делал узкое окно и сдвигал его относительно данных), то на 1/4 от начала частота 822.5897, на 1/2 - 822.5614, на 3/4 - 822.5633. При делении на 8 получается следующее:
822.5477
822.6066
822.5791
822.5578
822.5175
822.5213
822.5499
Видно, что чем мельче интервал анализа, тем больше разброс частоты. При большем усреднении разброс, естественно, уменьшается. Для данного файла болтанка частоты 9 единиц во втором знаке после запятой.
Я посмотрел на модели - синус плюс белый шум примерно равный тому, что есть в этом файле. Ошибка частоты - примерно 0.03 Гц максимум при делении на 8. Таким образом, здесь кроме шума действительно есть небольшое плавание частоты (примерно на те же 0.03 Гц). На полном интервале данный шум приводит к ошибке меньше 0.01 Гц. Отсюда следуют рекомендации (очевидные) - увеличивайте отношение сигнал/шум и боритесь с плаванием частоты.

[Ruslan1]

Видно, что чем мельче интервал анализа, тем больше разброс частоты. При большем усреднении разброс, естественно, уменьшается. Для данного файла болтанка частоты 9 единиц во втором знаке после запятой.
Я посмотрел на модели - синус плюс белый шум примерно равный тому, что есть в этом файле. Ошибка частоты - примерно 0.03 Гц максимум при делении на 8. Таким образом, здесь кроме шума действительно есть небольшое плавание частоты (примерно на те же 0.03 Гц). На полном интервале данный шум приводит к ошибке меньше 0.01 Гц. Отсюда следуют рекомендации (очевидные) - увеличивайте отношение сигнал/шум и боритесь с плаванием частоты.

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

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

[Alex11]

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

[Ruslan1]

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

на чистом синусе все шикарно, внутри бина ошибка (разность вычисленной и теоретической частот) соответствует предсказанной для гауссовского окна, именно так я и подбирал коэффициенты для окна (чтобы не более 0.001 Гц максимальная ошибка была).
Более того- и в приборе при подаче чистого синуса все отлично, ну прямо почти единицы наносекунд совпадают с частотомером.
А в реальности что-то не очень. Одни только очень серьезные палки на нечетных гармониках чего стоят. А фильтрануть не могу- полоса прибора позволяет куче гармоник пролазить. Реально хорошо вижу все нечетные гармоники, сколько аппаратный ФНЧ оставляет.

С шумами не смотрел на модели вообще, есть такой грех. Смотрел только разрешающую способность и точность при чистой синусоиде и радовался как слон Посмотрю... Кстати и с гармониками посмотрю....

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

Но это все так, рюшечки для себя любимого и набирание опыта с Матлабом. Поставленная при старте топика задача решена

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

[Alex11]

С гармониками бороться довольно просто - нужно, чтобы окно было выбрано так, чтобы на половине расстояния от основного тона до гармоники (по частоте) спектр спадал достаточно сильно. Если взять отношение е-7 - е-8, то достаточно. Хуже с шумом - он-то как раз в полосу попадает и всем мешает.
Производители суперкрутого прибора врут как могут. При соотношении сигнал/шум 21, как у них нарисовано, никакой точности в 3 знака после запятой получить нельзя. Какой там шум на самом деле - не видно, т.к. шкала должна быть логарифмическая.

[Ruslan1]

С гармониками бороться довольно просто - нужно, чтобы окно было выбрано так, чтобы на половине расстояния от основного тона до гармоники (по частоте) спектр спадал достаточно сильно. Если взять отношение е-7 - е-8, то достаточно. Хуже с шумом - он-то как раз в полосу попадает и всем мешает.
Производители суперкрутого прибора врут как могут. При соотношении сигнал/шум 21, как у них нарисовано, никакой точности в 3 знака после запятой получить нельзя. Какой там шум на самом деле - не видно, т.к. шкала должна быть логарифмическая.

Попробовал демпфировать датчик резистором- получил увеличение отношения сигнал/гармоника в 2 раза! Это я катушку датчика конкретно так нагружаю.... бум думать, тоже интересно...... Но не работает никто в токовом режиме, это я уже патент какой могу заявлять наверное
Гармоники проверил (на модели). Действительно зря боялся. На определение основной частоты при выбранном мной окне не влияют совсем: при 3-й,5-й 7-й гармониках с амплитудой равной основной гармонике рассчетная частота отличается от поданной примерно на 0.0006 Гц. То есть с этой стороны все в порядке.

А насчет той красивой картинки- так они точность 0.013% от значения обещают. даже для 100Гц это 0.013Гц, для более реальных 500Гц и выше- 0.065 Гц Не лучше чем у меня

[Alex65111]

Долго читал про потенциальную точность измерения, про то что здесь получается у народа (доли герца), но что-то с моими данными не могу сообразить.


1 вопрос. Имеются данные  dat1.zip ( 1.47 килобайт ) Кол-во скачиваний: 58
(после распаковки в матлабе загружаются load dat1). Как для данного конкретного случая (сигнал/шум, объем выборки, частота дискретизации 68кГц) оценить потенциальную точность измерения частоты? Формулу в предыдущих постах видел, но не пойму какое сигнал/шум подставлять, то ли во всей полосе, то ли только в полосе сигнала.

2. Насколько оптимально по точности измерение частоты следующим подходом, какую точность здесь можно ожидать?

Fs=68000;
N=512;
df=Fs/N;

sp=fft(s,N);

fmax=find(abs(sp(1:length(sp)/2))==max(abs(sp(1:length(sp)/2))));

f_ocen=(fmax(1)-1)*df;

f_s=f_ocen-df/1;
f_e=f_ocen+df/1;

m = 128*2;
w = exp(j*2*pi*((f_s-0)-(f_e+0))/(m*Fs));
a = exp(j*2*pi*(f_s-0)/Fs);

L1 = czt(s,m,w,a);
f_izm=find(L1==max(L1))*(2*df/m)+f_s;

[Ruslan1]

2. Насколько оптимально по точности измерение частоты следующим подходом, какую точность здесь можно ожидать?

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

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

[EvgenyNik]

Прежде чем экспериментировать и бросаться в пучину расчётов, рекомендую рассчитать два массива данных для интересующей длительности наблюдения и сравнить их.
1. Допустим, речь о частоте Fи
2. Желаемая точность dFи
3. Частота дискретизации Fд
4. Разрядность АЦП - N
5. Период наблюдения - T.
Таким образом, имеем T*Fд выборок, квантованных (в лучшем случае) по 2^(N-1) уровням...
Вычисляем 2 идеальных массива для частоты Fи и Fи + dFи и сравниваем их. В ряде случаев эти массивы (при чрезмерно оптимистичном подходе ) будут просто идентичными, тогда и копья ломать не стоит.
Если идти чуть дальше, то можно соотнести невязку этих двух массивов и одним из них с величиной, например, сигнал/шум реального АЦП. Это тоже покажет вам примерные идеальные возможности вашей системы.

[Dmitry Valento]

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

Ruslan1, подскажите пожалуйста адреса ресурсов, где вы нашли функции по вычислению окна Гаусса и гаусовской интерполяции.
К примеру, на сайте http://www.dsplib.ru/content/winadd/win.html указана немного другая формула по вычислению коэффициентов окна Гаусса.

Просто я попробовал применить ваш алгоритм на С++ на реальном сигнале от гитары, так он выдает всегда результат интерполяции на 0,5 больше чем исходная частота.
В математике я слаб, но ваш метод интерполяции мне показался уж очень простым

[Ruslan1]

Ruslan1, подскажите пожалуйста адреса ресурсов, где вы нашли функции по вычислению окна Гаусса и гаусовской интерполяции.

http://www.nicholson.com/rhn/dsp.html - там есть подборка ресурсов на тему уточнения FFT
http://www.mathworks.com/help/toolbox/signal/gausswin.html

К примеру, на сайте http://www.dsplib.ru/content/winadd/win.html указана немного другая формула по вычислению коэффициентов окна Гаусса.

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

Просто я попробовал применить ваш алгоритм на С++ на реальном сигнале от гитары, так он выдает всегда результат интерполяции на 0,5 больше чем исходная частота.
В математике я слаб, но ваш метод интерполяции мне показался уж очень простым

0.5 Гц? странно, о такой разнице я еще не слышал. Можно исходники и окна и интерполяции в студию? И файл с исходным оцифрованным сигналом, очень интересно.....
Если на 0.05 Гц- это уже ближе к телу....

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

[Alex11]

Правильное Гауссовское окно считается так:
sqrt(sqrt(PI)) / sqrt(sigma) * sqrt(FSIZE) * exp (-(x-FSIZE/2.0)*(x-FSIZE/2.0)/2.0/sigma/sigma);
PI - это число пи, FSIZE - размер Фурье в точках, sigma - параметр, задающий ширину окна.
При таком нормировочном коэффициенте накладывание окна не приводит к дополнительному коэффициенту в амплитудах.
По поводу гауссовской интерполяции - все правильно, но полезно вводить весовые коэффициенты, пропорциональные квадрату амплитуды сигнала - улучшает точность, т.к. шум при малом сигнале меньше влияет на определение положения вершины.

[Ruslan1]

Правильное Гауссовское окно считается так:
sqrt(sqrt(PI)) / sqrt(sigma) * sqrt(FSIZE) * exp (-(x-FSIZE/2.0)*(x-FSIZE/2.0)/2.0/sigma/sigma);
PI - это число пи, FSIZE - размер Фурье в точках, sigma - параметр, задающий ширину окна.
При таком нормировочном коэффициенте накладывание окна не приводит к дополнительному коэффициенту в амплитудах.
По поводу гауссовской интерполяции - все правильно, но полезно вводить весовые коэффициенты, пропорциональные квадрату амплитуды сигнала - улучшает точность, т.к. шум при малом сигнале меньше влияет на определение положения вершины.

Извиняюсь за очередной глупый вопрос: а как выбрать ширину окна (sigma)?
Например, хочу задавить третью гармонику, которая в худшем случае находится на расстоянии 105 бинов, если рассматривать частотную область. Это значит, что мне нужно взять окно, которое в временной области задавит находящийся на расстоянии 105 семплов от ценрального отчета (FSIZE/2) ? Или все совсем не так?
У меня получилось что при сигма=43 имею подавление 60 dB во временной области.

Или достаточно, чтобы к границам окна наблюдения (отчеты 0 и FSIZE) упало до нуля ( ну или -60 dB хотя бы) ?

[Alex11]

Ширина выбирается в зависимости от задачи. Мне нужно было, чтобы вклад второй гармоники в амплитуду основного сигнала не превышал 1е-8, ну я и подбирал сигму так, чтобы посередине между основным тоном и второй гармоникой был провал не хуже 1е-8. Для определения только частот можно, думаю, ограничиться разделением на уровне 1е-6.
По поводу соотношения окна в частотной и временной областях - что-то там похожее должно быть, но я не исследовал точно. Может какой-то коэффициент вылезти. Лучше смотреть сразу в спектре.

[Dmitry Valento]

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

[Ruslan1]

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

да пожалуйста.

Ширина выбирается в зависимости от задачи. Мне нужно было, чтобы вклад второй гармоники в амплитуду основного сигнала не превышал 1е-8, ну я и подбирал сигму так, чтобы посередине между основным тоном и второй гармоникой был провал не хуже 1е-8. Для определения только частот можно, думаю, ограничиться разделением на уровне 1е-6.
По поводу соотношения окна в частотной и временной областях - что-то там похожее должно быть, но я не исследовал точно. Может какой-то коэффициент вылезти. Лучше смотреть сразу в спектре.

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

[Alex11]

Шкалу вертикальную сделайте линейную, а не логарифмическую - тоже ничего не увидите. У них же не децибелы, а амплитуда.

[Alexey Lukin]

sqrt(sqrt(PI)) / sqrt(sigma) * sqrt(FSIZE) * exp (-(x-FSIZE/2.0)*(x-FSIZE/2.0)/2.0/sigma/sigma);
При таком нормировочном коэффициенте накладывание окна не приводит к дополнительному коэффициенту в амплитудах.

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

[Alex11]

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

[Fetronics]

Вопросы о разрешающей способности измерения частоты за доли секунды с разрешением меньше Герца можно посмотреть здесь, очень невероятные вещи узнаете!
http://electronix.ru/forum/index.php?showt...mp;#entry703059

[rudy_b]

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

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

1. Частота колебаний зависит от их амплитуды. Грубо говоря, если струна колеблется с большей амплитудой - она больше растянута и ее жесткость выше. Единственный выход - минимизировать амплитуду ((или стабилизировать или выбирать участок с определенным ее значением).

2. В струне, особенно не очень корректно закрепленной, присутствуют колебания на нескольких очень близких частотах. Это связано, например с тем, что струна имеет чуть разные длины и жесткости для вертикального и горизонтального колебания, например на гитаре. А горизонтальная и вертикальная моды переходят одна в другую, в результате появляется медленная модуляция амплитуды на фоне общего затухания сигнала.

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

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

[ys05]

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

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

[rudy_b]

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

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

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

[tmtlib]

Если кому ещё интересна эта тема, вот хороший материальчик:
http://stackoverflow.com/questions/4633203...-between-frames

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

Код

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

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

[fontp]

Если кому ещё интересна эта тема, вот хороший материальчик:
http://stackoverflow.com/questions/4633203...-between-frames

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

CODE

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

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

Хороший алгоритм по типу штангенциркуля. Только Вы забыли сказать, что такое k. Из-за этого, если не читать матерьяльчик, то непонятно, что здесь написано. Здесь k - это номер бина в котором максимум энергии. В матерьяльчике автор показывает, что правильную частоту можно получить не только по этому "максимальному" бину но и по FFT_N/FFT_STEP соседних. Поэтому не возникает проблем с частотами на границе между бинами при перекрытии блоков FFT. Без нахлёста с такими пограничными частотами алгоритм работать не будет, шумы не дадут. Систематическая ошибка у алгоритма отсутствует. Интересно посмотреть как алгоритм ведёт себя по отношению к шуму - судя по всему выйдет на точность оценки максимального правдоподобия CRLB, особенно если ещё и научиться делать взвешенную оценку по нескольким значениям возле максимума (остальные FFT_N/FFT_STEP оценок из-за малой энергии бинов реально утонут под шумом)

[ivan219]

Если кому ещё интересна эта тема, вот хороший материальчик:
http://stackoverflow.com/questions/4633203...-between-frames

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

Код

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

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

А можно этот код расписать иначе? Так как в матлабе не силён. Если можно то блок схемой или просто словами.

[tmtlib]

А можно этот код расписать иначе? Так как в матлабе не силён. Если можно то блок схемой или просто словами.

Я планирую выложить примерчик с интересными методами здесь. Но думаю это будет не скоро. А пока попробую на словах, например, для фурье на 1024 отсчёта:
1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].
2) Для каждого из 512 элемента считаем фазу. Если по простому, то фаза[i]:=ArcTan(Im[i] / Re[i]); Можно сделать получше, расписав разные случаи (если Re[i]=0 и т.д.)
3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[i], т.е. фаза[i] от текущего FFT. Заметил, что метод лучше работает если применить сглаживающее окно.

Забыл написать: например, константа SAMPLE_RATE=44100, FFT_N =1024, FFT_STEP = 100 (количество точек, на которые смещаем исходные данные. 100 намного меньше 1024, поэтому результат лучше)

Сообщение отредактировал tmtlib - Nov 13 2011, 14:06

[ivan219]

И всё таки не много не понятно.
Что делает эта функция?

Код

delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval

[thermit]

вычисляет остаток от деления дельты на два пи, вероятно...

Сообщение отредактировал thermit - Nov 13 2011, 21:17

[tmtlib]

И всё таки не много не понятно.
Что делает эта функция?

Код

delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval

Фазу приводят к диапазону -3.14...+3.14, пример в градусах: это что-то наподобие того, что 350 градусов привратится в -10 и т.д. Т.е. вектор (Im,Re) всегда в правом полукруге от -90 до 90 градусов.

[fontp]

Фазу приводят к диапазону -3.14...+3.14, пример в градусах: это что-то наподобие того, что 350 градусов привратится в -10 и т.д. Т.е. вектор (Im,Re) всегда в правом полукруге от -90 до 90 градусов.

Именно на этом этапе алгоритм перестанет работать, если нет перекрытия. Как показывает автор "матерьяльчика" в последней табличке (Pass 4) для отсутствия перекрытия.
http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/
Для частот на границе бинов получаем 2 значения частоты для примерно одинаковых величин спектральной мощности, и только одно из них правильное. Поскольку максимум двух равных величин при реальных вычислениях определяется шумом, то в этом случае алгоритм будет выбирать правильное значение с вероятностью 0.5. Впрочем можно попытаться определить частоту другим методом по каждому из блоков отдельно и потом взять то значение которое ближе к этой оценке, но это уже как бы гибрид

2) Для каждого из 512 элемента считаем фазу. Если по простому, то фаза[i]:=ArcTan(Im[i] / Re[i]); Можно сделать получше, расписав разные случаи (если Re[i]=0 и т.д.)
3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[i], т.е. фаза[i] от текущего FFT. Заметил, что метод лучше работает если применить сглаживающее окно.

Для определения частоты вам не нужны все FFT_N значений фазы, поскольку только FFT_N/FFT_STEP значений фазы около максимума спектральной мощности правильные, а остальные искажены врэпингом на 2*pi (Pass 5). Более того в реальной ситуации присутствуют шумы и большинство из этих FFT_N/FFT_STEP правильных оценок (если их много) будут реально сильно зашумлены из-за очень быстрого падения уровня сигнала в бине при отклонении от центрального. Мы помним, что отклик бинов на синусоиду падает в соответствии с функцией спектрального окна, и в лучшем случае он будет спадать как синк или быстрее. В реальной ситуации только 2-4 центральных бина будут давать значимую оценку частоты, причем наименее искажены шумом будут только одна (если частота близко к центру бина) или две (если частота на границе бинов) оценки для центральных бинов (где максимум спектральной мощности). Остальные "утонут" под шумом
Все 512 значений могут представлять интерес только для визуализации "сонограм", поскольку они, в большинстве, шумоподобны.
Поэтому приведеный алгоритм "недоделаный" - мало того что он неэффективен, но он не содержит последнего шага, поскольку он генерирует FFT_N значений частоты из которых в лучшем случае (отсутствие шума) только FFT_N/FFT_STEP правильны и не сказано, как правильные выбрать из массива оценок, большинство из которых неправильны, см. таблички у автора
http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/
По смыслу понятно, что выбирать нужно где-то возле максимума спектральной мощности

С учетом сказанного простейший алгоритм будет выглядеть так:
1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].
2) Вычисляем энергию для каждого k E[k]=Re*Re+Im*Im
3) Находим максимум E[k] - пусть будет E() для kmax. Здесь вообще-то нужно проверять что kmax не прыгает. Если изменение kmax от фрейма к фрейму больше 1 по модулю - то это не синусоида, выход по ошибке
4) Для элементов kmax, kmax-1, kmax+1 считаем фазу. фаза[]:=ATAN2(Im[], Re[]); другие элементы на самом деле не нужны, в них отсутствует полезная информация
3) храним посчитаные фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[kmax], т.е. фаза[kmax] от текущего FFT. Поскольку kmax гарантировано не меняется больше чем на 1, то в любом случае значение фазы для kmax посчитаны и для предыдущего фрейма и есть в массиве

Если перекрытие FFT_N/FFT_STEP >2 то оценку можно немного улучшить если использовать не одно значение максимума, а два первых максимума и построить взвешенную оценку для частоты по этим двум оценкам.
Если главный максимальный бин k0=kmax, а следующий по энергии k1=arg max(E(kmax-1), E(kmax+1)), то можно построить оценку типа
(f(k1)*E(k1)+f(k0)*E(k0))/(E(k0)+E(k1)) с большим иммунитетом к шуму, главным образом для частот на границе бинов (поскольку в этом случае будет две равноценных оценки)

А так да, идея алгоритма на вскидку смотрится очень хорошо

[alex_os]

...
(f(k1)*E(k1)+f(k0)*E(k0))/(E(k0)+E(k1)) с большим иммунитетом к шуму, главным образом для частот на границе бинов (поскольку в этом случае будет две равноценных оценки)

f(k1) это измеренная частота для бина k1 ? А что если приращение фазы измерять по нескольким бинам в окрестности
максимума? Взвешивание получится автоматически.

Код

        delta =  angle( X_1(kmax-n:kmax+n)' * X(kmax-n:kmax+n) )

X_1 - предидущий выход FFT (вектор столбец),
X - текущий выход FFT,
kmax - индекс бина с максимальной энергией,
n - число 0...3

[fontp]

f(k1) это измеренная частота для бина k1 ? А что если приращение фазы измерять по нескольким бинам в окрестности
максимума? Взвешивание получится автоматически.

CODE

        delta =  angle( X_1(kmax-n:kmax+n)' * X(kmax-n:kmax+n) )

X_1 - предидущий выход FFT (вектор столбец),
X - текущий выход FFT,
kmax - индекс бина с максимальной энергией,
n - число 0...3

Можно бы и по всем. только не все бины одинаково полезны. Поэтому это будет не оптимально.
Усреднять надо в соответствии с уровнем сигнала каждого бина. Вычисление фазы в бинах с низким уровнем сигнала по отношению к шуму бесмысленны - это уже фаза компоненты шума а не сигнала. А уровень сигнала убывает при отклонении от центральных бинов достаточно быстро (быстрее чем 1/n), в то время как шум остается постоянным. Поэтому при низком уровне сигнал/шум лучше усреднять по 2-3-м, причем с учетом энергии (можно даже расписать точные формулы, выписывая плотность вероятности ошибки фазы в бине как функцию отношения сигнал/шум - это можно найти в учебниках по спектральному анализу или статрадиофизике - и максимизмизируя произведение вероятностей).

При очень высоком уровне сигнал/шум, хотя бы 20дб, действительно усреднять вроде можно по всем четырём (FFT_N/FFT_STEP в общем случае перекрытия, но тоже с учетом энергии бина), поскольку ошибки измерения фазы статистически независимы, но различны для разных бинов, в зависимости от энергии

да, f(k) - измеренная частота для бина k. Я не максимизировал функцию вероятности, а просто написал на вскидку формулу, которая будет вести себя оптимально для частот и посередине бина (f(k0)) и на границе бинов (f(k0)+f(k1))/2,
если этого не сделать, то точность упадёт на границе

ЗЫ. Вообще-то вычисление фазы через квадратуры всегда усиливает шум и помехи. Поэтому при вычислении частоты через фазу не стоит пренебрегать использованием спектральных окон, как в примере - чтобы подавить хвосты от возможных синусовых помех. Они и всегда мешают, а этому методу будут особенно

[alex_os]

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

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

[fontp]

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

Виноват. Про angle (в Матлабе?) не знал, думал среднее арифметическое фазы.

[alex_os]

Виноват. Про angle (в Матлабе?) не знал, думал среднее арифметическое фазы.

Сорри, не написал что матлаб. angle -аргумент комплексного числа.

[rudy_b]

Можно бы и по всем. только не все бины одинаково полезны...

Более правильно - по МНК подбирается параметры гаусса (огибающей) (при соответствующей весовой функции) с весами, пропорциональными мощности. При этом точность возрастает на порядок.

[fontp]

Более правильно - по МНК подбирается параметры гаусса (огибающей) (при соответствующей весовой функции) с весами, пропорциональными мощности. При этом точность возрастает на порядок.

Так предложили уже взвешивать, если говорить о алгоритме с разностью фаз ДПФ перекрывающихся фреймов.
Или Вы о чём то своём?

Сообщение отредактировал fontp - Nov 16 2011, 11:10