bahurin

Ну как бы они все стянуты в ноль.

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

Согласен при равной разрядности аккумулятора интегратора и дифференциатора все будет хорошо

Допустим есть каскад I-C, где I интегратор, а С дифференциатор H(z) =1-z^-1. Допустим аккумулятор интегратора 4 бита. На вход поступают одни единицы 1 1 1 1 1 1... По идее такой каскад на выходе ничего не должен менять, потому что числитель H(z) равен знаменателю, но при 4 битной разрядности интегратора когда он досчитает до 15 на следующем такте он переполнится и выдаст 0. Дифференциатор выдаст разницу 0-15 = - 15 и потом снова 1 1 1 1... и через 15 тактов снова на выходе - 15. Если вы используете дециматор, то он отбросит эту левую -15, а если каскад просто используется как фильтр скользящего среднего, то в выходном сигнале - 15 пойдёт дальше как есть.

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

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

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

Cic filter имеет конечную импульсную характеристику

Если без высшей математики, то рекурсивных фильтров с линейной фчх не бывает. Можно лишь аппроксимировать линейную фчх постановкой allpass корректоров.

Думаю что проблема в использовании IFFT. Дело в том, что для IFFT нужны равноотстоящие спектральные отсчеты в диапазоне от 0 до Fs/2 Гц. Скорее всего в вашем файле они от некоторой частоты F0 до F1. В этом случае надо дополнить нули в массив чтобы он стал от 0 до Fs/2 Гц, потом сделать из него массив симметричный относительно Fs/2 (чтобы результат IFFT стал реальным). И только после этого надо применять IFFT. но даже в этом случае не уверен, что будет правильный результат. Надо смотреть данные с вашего нетворк анализатора. Возможно надо применять не IFFT, а напрямую численное интегрирование обратного преобразования Фурье.

Что значит из файла с s параметрами? Откуда файлы?

После ifft при обнуленной ФЧХ обязательно надо делать перестановку импульсной характеристики. Нельзя прямо сразу сверку с результатом ifft! Ещё желательно окошком после перестановки.

Не понял почему вы считаете приведенный пример неустойчивым? Вывел на печать параметры g1 и g2: g1 = 0.0319 g2 = 0.0010. Все в пределах зоны устойчивости. также не совсем ясно что вы подразумеваете под границами применимости конкретной данной ФАПЧ? Каждая следящая система должна проектироваться под конкретные параметры. Не всегда можно взять готовое решение и прикладывать для своей задачи. Так если тактовую поменяете, надо пересчитать петлевой фильтр (он должен стать более узкополосным чтобы соблюсти полосу захвата).

Эта библиотека все умеет и fir и iir. Документация тут Я там правда не все задокументировал, но если что, пишите помогу.

Можно ещё посоветовать посмотреть в сторону python. Octave практически повторяет синтаксис матлаба. Надо только нужные тулбоксы включать например load pkg signal

a[0] = z[2]; a[3] = DSPL_FARROW_LAGRANGE_COEFF*(z[3] -z[0]) + 0.5*(z[1] - z[2]); a[1] = 0.5*(z[3] - z[1])-a[3]; a[2] = z[3] - z[2] -a[3]-a[1]; res = polyval(a, 3, &x, 1, (*y)+k); Если z обозначить как y то отличий ровно 0. Только приведённые обертки позволяют ещё сдвигать на дробное время и ресаэмплировать в произвольное отношение P/Q. Кроме того есть ещё реализация на кубических сплайнах, которая не требует коэффициента 1/6

Если кому интересно есть Сишная реализация фарроу фильтров

Вот наглядный пример почему так делать не стоит. :laughing:

А почему при 16 битном представлении у вас максимальный к-т фильтра 2577? Что мешает усилием воли туда 32767 записать? Надо понять где точка в ваших к-тах и столько разрядов по выходу фильтра сделать раунд для корректного масштабирования

Комплексная экспонента есть exp(-jwt) = cos(wt)-j*sin(wt). Фаза есть atan(imag(S)/real(S)) . Соответсвенно для косинуса получим фазу 0 или pi поскольку спектр будет чисто вещественным, а для синуса спектр будет чисто комплексным и фаза будет +-pi/2

если резулатат вашего алгоритма Герцеля отличается от fft то значит где-то ошибка. Результат должен совпадать. Теперь по существу. Вот матлаб код для проверки. Берем 500 отсчетов синусоиды s = sin(2*pi*50*t/Fs) на частоте 50 Гц и рассчитываем 50-ый отсчет FFT используя алгоритм герцеля. Получаем фазу -90 градусов, поскольку фаза 0 будет если взять s = сos(2*pi*50*t/Fs). Если посмотрите спектр FFT, то нет эффекта растекания поскольку я взял 500 отсчетов. Если вы берете 512 отсчетов то эффект растекания будет и фаза изменится. При этом если посмотрите из примера разница между отсчетом FFT и герцелем ничтожно мала (1E-12). Так и должно быть поскольку алгоритм Герцеля и рассчитывает один отсчет FFT. clear all; close all; clc; %% input data N = 500; t = 0:N-1; Fs = 500; k = 50; s = sin(2*pi*50*t/Fs); %% Goertzel algorithm W = exp(2i*pi*k/N); v = zeros(1,N); a = 2*cos(2*pi*k/N); v(1) = s(1); v(2) = s(2) + a*v(1); for n = 3:N v(n) = s(n) + a*v(n-1) - v(n-2); end G = W*v(n) - v(n-1); %% fft S = fft(s); figure; plot(abs(S)) %% fprintf('разница между fft и Герцелем: %e\n', G - S(k+1)); fprintf('фаза: %.4f\n', angle(G)*180/pi); PS по поводу накопления ошибок поставьте в исходных данных %% input data N = 500000; t = 0:N-1; Fs = 500; k = 50000; и все равно ошибка будет на уровне 1E-6 при абсолютном значении спектра 1E6 (т.е. на 12 порядков ниже)

Почему нет соотвествия по периоду? Частота сигнала 1 Гц палка на частоте 1 что не так? Амплитуда просто должна быть поделена на N если хотите получить честное значение 0.5. В скрипте 6 строчек с какой именно вы не согласны?

Открываем матлаб вставляем и запускаем скрипт: N = 64; t = (0:N-1)/N; s = sin(2*pi*t); subplot(211), stem(t,s), grid, title('один период синусоиды'); S = abs(fft(s)); subplot(212), stem(0:N-1,S), grid, title('спектр одного периода синусоиды'); смотрим график: Очевидно что на графике только одна палка S(1) все остальные равны нулю (надеюсь не надо объяснять что такое вторая палка). После этого извиняемся перестаем дерзить и идем допиливать свою программу до рабочего состояния или перестаем называть свою программу анализатором спектра поскольку к спектральному анализу она в своем текущем виде не имеет никакого отношения.

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

1. Почему не увенчалось успехом? Не удалось сделать ifft в матлабе и умножить на окно? 2. Увеличивая длину фильтра и применяя окно Кайзера вы сможете добиться требуемых характеристик.