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

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

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

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

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

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



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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Оставлю на вашей совести, обвинение Марпла в публикации неработающего метода в книге. Говорю же в книге, особенно переводной, могут быть опечатки, особенно типично I пропечатать как 1 или J, могут быть пропущены даже строки в программе.

Ладно, помогу Вам немного, раз уж ввязался. В книге есть ссылки на публикации по каждому методу, и в большинстве публикаций в то время было принято приводить псевдокод и, часто, код на ФОРТРАНе. А ковариационный метод для Марпла "проприетарный", в том смысле, что это метод Марпла, как метод Левинсон-Дурбина - это метод Левинсона и Дурбина, а не Берга. Поэтому он особенно постарался.
Ковариационный метод Марпл предлагал даже в более общем и интересном варианте -
для идентификации FIR-систем по короткой последовательности (типа статический эквалайзер модема по известной преамбуле). То есть Yn = FIR(Xn-1), а не Xn=FIR(Xn-1) как в AR
Эта работающая реализация не настолько стара, чтобы быть мне недоступной.
Вот эта статья c программой (FORTRAN в самом крутом, комплексном варианте), сравните со своей реализацией, заменяя вектор Y на сдвинутый Х, получится ваша AR-модель.
ЗЫ.В принципе у меня есть РАБОТАЮЩАЯ реализация на С, но не очень хочется заморачиваться с ее размещением (С-файл в форуме не атачится), да Вы и не просили, да и вообще, если программа не понравится - скажите какую-то гадость - лучше не буду навязываться

Спасибо большое за помощь!

Можно сделать дискретное преобразование Фурье (DFT) на неравномерной сетке. Оно немного (зависит от степени неравномерности) хуже обусловлено, и требует большего объёма вычислнеий по сравнению с DFT (нужно решать систему линейных уравнений)



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

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

Я бы частоту самплинга повыше выбрал. Вы ставите anti aliasing filter? Он частоту Найквиста (Котельникова) подавит, что я думаю вам нежелательно. Не ставить фильтр также нежелательно. Без фильтра можно получить "чудеса". Продолжительность измерения определяет разрешение шкалы результата. Допустим вы измеряли 1000 секунд, тогда после преобразования Фурье разрешение будет 1/1000 = 0.001 Герц. Конечно можно и получить ниже, но выше не получится.




Для одной частоты не надо делать преобразование Фурье. Достаточно считать корелляции с синусом и косинусом заранее известной частоты.

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

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



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

Все определяется целью, в АОНАх или DTMF детекторах этот алгоритм блестяще справляется.

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




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

Действительно не даст из-за систематического смещения оценки.Но

1. Это смещение очень мало для зашумленного сигнала ( хуже 30-40 дб) по сравнению с случайной ошибкой, обусловленной шумом. Причем эта случайная ошибка достигает минимальных предельных значений, достижимых любым другим способом, поскольку она приближается к оценке максимального правдоподобия, которая известна как предельная оценка Крамера-Рао

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

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

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