Спектральный анализ на сверхнизких частотах Опции

[Crowbar]

Допустим, требуется получить разложение спектра частот до 20Гц с точностью 0,01Гц и выше. Каким образом это реализуется, помимо самого простого способа, как поставить частоту отцифровки на 40Гц, выставить кол-во сэмплов на 4000 и ждать больше полутора минут завершения очередного цикла?

[mikalaha]

Допустим, требуется получить разложение спектра частот до 20Гц с точностью 0,01Гц и выше. Каким образом это реализуется, помимо самого простого способа, как поставить частоту отцифровки на 40Гц, выставить кол-во сэмплов на 4000 и ждать больше полутора минут завершения очередного цикла?

Не обязательно. Просто набираете некоторое количество отсчетов (например 80 - 2 секунды записи). Затем производите БПФ требуемой размерности (например, 4096 - точность как раз около 0.01Гц). Недостающие до 4096 отсчеты (4016) заполняете нулями. И выполняете БПФ.

[Alex255]

Не обязательно. Просто набираете некоторое количество отсчетов (например 80 - 2 секунды записи). Затем производите БПФ требуемой размерности (например, 4096 - точность как раз около 0.01Гц). Недостающие до 4096 отсчеты (4016) заполняете нулями. И выполняете БПФ.

То есть Вы хотите получить информацию о процессе периодом 100сек за 2секунды? Увы...

To Crowbar: Да, именно так как Вы сказали

[mikalaha]

То есть Вы хотите получить информацию о процессе периодом 100сек за 2секунды? Увы...

To Crowbar: Да, именно так как Вы сказали

На самом деле все отлично работает.
Например:
сигнал - тон с частотой 10 Гц (период 0.1 сек) - в одной секунде 10 периодов!!
при записи с дискретизацией 40 Гц одной секунды получаем 40 отсчетов, в которых содержится
целых 10 периодов тона.
Если теперь дополнить нулями и взять FFT мы получим этот тон.

Прикреплены картинки
time - сигнал во временной области;
freq - после FFT 4096
freq2 - увеличенная частотная позиция 10 Гц после FFT (что и требовалось)

На самом деле все отлично работает.
Например:
сигнал - тон с частотой 10 Гц (период 0.1 сек) - в одной секунде 10 периодов!!
при записи с дискретизацией 40 Гц одной секунды получаем 40 отсчетов, в которых содержится
целых 10 периодов тона.
Если теперь дополнить нулями и взять FFT мы получим этот тон.

Прикреплены картинки
time - сигнал во временной области;
freq - после FFT 4096
freq2 - увеличенная частотная позиция 10 Гц после FFT (что и требовалось)

То же самое наблюдается и при тоне с частотой например (10.05 Гц). Единственное различие вышло следующее: без прореживания нулями пик после Фурье вышел на частоте 10.0488 Гц, а с прореживанием - 10.0399 Гц

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

 

[fontp]

Дополнение нулями работает, но только для одиночного тона. На самом деле вы получаете грубую сетку частот для оценки по БПФ, а более тонкие частоты получили интерполяцией (дополнение нулями означает некоторую интерполяцию промежуточных точек БПФ). Для оценки произвольного спектра это работать, очевидно, не будет. Попробуйте, например, разрешить два тона, отстоящих друг от друга на 0.01 гц - 10гц и 10.01 гц. Разрешить означает, что они присутствуют в спектре одновременно.
Неуспех гарантирован принципом неопределённости

[el34]

2 Crowbar возмите - в инете есть книжка Марпла


С.Л. Марпл-младший. Цифровой спектральный анализ и его приложения. (гл. 1-7)
С.Л. Марпл-младший. Цифровой спектральный анализ и его приложения. (гл. 8-16)
http://dsp-book.narod.ru/books.html

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

[Crowbar]

О, а я ее чего-то пропустил, когда там был, спасибо.

[mikalaha]

Дополнение нулями работает, но только для одиночного тона. На самом деле вы получаете грубую сетку частот для оценки по БПФ, а более тонкие частоты получили интерполяцией (дополнение нулями означает некоторую интерполяцию промежуточных точек БПФ). Для оценки произвольного спектра это работать, очевидно, не будет. Попробуйте, например, разрешить два тона, отстоящих друг от друга на 0.01 гц - 10гц и 10.01 гц. Разрешить означает, что они присутствуют в спектре одновременно.
Неуспех гарантирован принципом неопределённости

Согласен, погорячился. Надо подумать об обучении

[Waso]

Занимаюсь получением спектра по модифицированному ковариационному алгоритму, описанному в книге Марпла мл., на которую тут ссылались.
Написал прогу. Вернее списал с книжки, перевел с тамошнего фортрана на экселевский бейсик (ну надо оно мне там!!).
Но она считает неправильно. В конце той книги должна быть тест-последовательность, по которой можно отладить все процедурки,
приведенные в книге. Но в моем дэжавюшнике книга обрывается на стр 547, где еще идет окончание 16-й главы... Откуда я только не качал
эту книженцию - везде лежит одно и тоже. Иногда по половинкам иногда целая, но косяки все теже.
У кого есть ПРИЛОЖЕНИЕ II с 64-точечной комплексной тест-последовательностью марпла - помогите пожалуйста! Можно прямо тут положить.
====================================
А все. Сам нашел. )))) Там-же есть постраничная выкладка этой книги. Там все от начала и до конца и в хорошем качестве. Только чтоб показывало - надо скачать 0584.djbz обязательно. и можно по страничкам смотреть.

[ne_ya]

Занимаюсь получением спектра по модифицированному ковариационному алгоритму, описанному в книге Марпла мл., на которую тут ссылались.
Написал прогу. Вернее списал с книжки, перевел с тамошнего фортрана на экселевский бейсик (ну надо оно мне там!!).
Но она считает неправильно. В конце той книги должна быть тест-последовательность, по которой можно отладить все процедурки,
приведенные в книге. Но в моем дэжавюшнике книга обрывается на стр 547, где еще идет окончание 16-й главы... Откуда я только не качал
эту книженцию - везде лежит одно и тоже. Иногда по половинкам иногда целая, но косяки все теже.
У кого есть ПРИЛОЖЕНИЕ II с 64-точечной комплексной тест-последовательностью марпла - помогите пожалуйста! Можно прямо тут положить.
====================================
А все. Сам нашел. )))) Там-же есть постраничная выкладка этой книги. Там все от начала и до конца и в хорошем качестве. Только чтоб показывало - надо скачать 0584.djbz обязательно. и можно по страничкам смотреть.

У меня та же самая проблема -- не подскажите, как решили? До того реализовывал метод Берга -- рекурсивное нахождение дисперсии белого шума (rho) давало неправильные результаты. После того, как рекурсию заменил на "грубую силу", начал получать правильный спектр. Может, здесь тоже с этим проблема? Или я где-то с коэффициентами напутал?

[fontp]

У меня та же самая проблема -- не подскажите, как решили? До того реализовывал метод Берга -- рекурсивное нахождение дисперсии белого шума (rho) давало неправильные результаты. После того, как рекурсию заменил на "грубую силу", начал получать правильный спектр. Может, здесь тоже с этим проблема? Или я где-то с коэффициентами напутал?

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

Если данных мало,среди блочных методов ковариационные методы лучше всего.

Скользящая же рекурсия (алгоритм Fast RLS) работает условно - то есть он сходится до определенного предела, после чего обязательно разваливается, поскольку вблизи точного решения обычно алгоритм становится сингулярным, выходит на границу устойчивости.
Самый простой способ, используемый для борьбы с этим - это периодическая реинициализация, там сказано. Кроме того существуют более устойчивые варианты FRLS, чем классический, приведенный у Марпла (как адапитивная Калмановская AR-модель)

типа этого FRLS, устойчивого в большинстве случаев
http://www.wcl.ece.upatras.gr/CSNDSP//cont...0Audio/A8.1.pdf

[ne_ya]

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

Если данных мало,среди блочных методов ковариационные методы лучше всего.

Спасибо за ответ.

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

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

[fontp]

Спасибо за ответ.

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

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

Я когда-то реализовывал все блочные методы, но очень давно. Когда еще был ФОРТРАН. Деталей не помню. Там в тексте программ главное не путать I и 1, и не запутаться в индексах, если перекладывать на другой язык

[_pv]

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

1) если просто посчитать интеграл Фурье, какая получится полоса у полученного спектрального отсчёта?
2) как можно простым образом эту полосу контролирумо уменьшить?

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

[ne_ya]

Я когда-то реализовывал все блочные методы, но очень давно. Когда еще был ФОРТРАН. Деталей не помню. Там в тексте программ главное не путать I и 1, и не запутаться в индексах, если перекладывать на другой язык

откомпелировал код на фортране, запустил -- ответ с предоставленным в книге не совпадает. Зато очень даже похож на то, что при переписывании на c# получалось. Так что все-таки вряд ли "все блочные методы у Марпла рабочие", к сожалению.