Оставлю на вашей совести, обвинение Марпла в публикации неработающего метода в книге. Говорю же в книге, особенно переводной, могут быть опечатки, особенно типично 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- линейной системы по минимуму квадратов быстрым алгоритмом (импульсного отклика, но не спектрального)