[stealthisname 7]

9 hours ago, C2000 said:

Еще не подскажите как в Матлабе вот эти вильтры получить:

Все пропускающие фильтры

 

В fdatool таких не нашёл

в fdatool есть генерация всепропускающих фильтров для формирования необходимой групповой задержки

и эти фильтры формируются из цепочкой фильтров, аналогичных фильтрам  рассматриваемым в указанной статье

применительно к задаче по формированию определенного ФЧХ использовать напрямую эту функцию fdatool не получится

 

так как в статье указаны результаты расчеты передаточных характеристик, можно получить именно эти фильтры

например для начала статьи, для фильтров первого порядка

для одного фильтра с заданной k

k = -.9;

b = [k,1];
a = [1,k];

hfvt = fvtool(b,a);
hfvt.Analysis = 'phase';
hfvt.PhaseUnits = 'degrees';

 

для серии фильтров

k_arr = -.9:.1:.9;

c = {};
for i = 1:numel(k_arr)
    k = k_arr(i);
    b = [k,1];
    a = [1,k];
    c = [c,b,a];
end

hfvt = fvtool(c{:});
hfvt.Analysis = 'phase';
hfvt.PhaseUnits = 'degrees';

 

[C2000 2]

Нашел еще вот такую интересную статью:

hilbert_Olli_Niemitalo (90 degree phse shift filter).pdf

Вот только автор не рассказывает как он коэффициенты к этим фильтрам подбирал, только анализ результата.

Есть ли какие то мысли, как эти фильтры рассчитывались?

 

Там в конце дана формула одного звена:  out(t) = a^2*(in(t) + out(t-2)) - in(t-2)

Как получить АЧХ и ФЧХ именно по функции, в конечном итоге хотелось бы Си файл с функцией фильтра подцеплять и проверять что получается.

[C2000 2]

 

....

[C2000 2]

22 hours ago, stealthisname said:

 

% комплексные частотные харатеристики
h  = sum(b .'.*exp(-1i*2*pi*w.*(1:N+1).'));
h0 = sum(b0.'.*exp(-1i*2*pi*w.*(1:N+1).'));

 

Можно здесь подробнее, что означают эти иероглифы)))

Хотел то же самое проделать с IIR фильтром, но как, не понял пока что.

[stealthisname 7]

2 hours ago, C2000 said:

Можно здесь подробнее, что означают эти иероглифы)))

Хотел то же самое проделать с IIR фильтром, но как, не понял пока что.

эти вычисления -  подстановка в передаточную функцию H(z) частоты z = exp(j*2*pi*w)

для БИХ фильтров вычисления аналогичные

с оформлением в виде функции и примером использования:

% КЧХ для КИХ
func_h_fir = @(b,w) sum(b.'.*exp(-1i*2*pi*w.*(0:numel(b)-1).'),1);

% КЧХ для БИХ
func_h_iir = @(b,a,w) func_h_fir(b,w)./func_h_fir(a,w);

% АЧХ, дб
func_h_mg = @(b,a,w) mag2db(abs(func_h_iir(b,a,w)));

% ФЧХ, °
func_h_ph = @(b,a,w) rad2deg(unwrap(angle(func_h_iir(b,a,w))));

% цифровая нормированная частота
w = 0:1e-3:1;

% коэффициенты передаточной функции
b1 = [1,1]/2;
a1 = [1];

b2 = [0,.1];
a2 = [1,-(1-.1)];

% частотные характеристики
h1_mg = func_h_mg(b1,a1,w);
h1_ph = func_h_ph(b1,a1,w);

h2_mg = func_h_mg(b2,a2,w);
h2_ph = func_h_ph(b2,a2,w);

% график АЧХ
plot(w,h1_mg,w,h2_mg);

% график ФЧХ
plot(w,h1_ph,w,h2_ph);

 

[C2000 2]

Подскажите как в Matlab посчитать вот это:

 

При использовании symsum с аргументов inf Матлаб виснет.

[stealthisname 7]

54 minutes ago, C2000 said:

Подскажите как в Matlab посчитать вот это:

 

 

один из возможных способов расчета подобных выражений рассмотрим на примере суммы в числителе выражения

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

рассмотрим расчет для первой сотни слогаемых

q = .99;
m = 0:100;
cs = cumsum( q.^(m.*(m+1)));
plot(m,cs);
grid on;

построим для разных значений q

q = [0:.1:.9,.95,.99];
m = 0:100;
cs = cumsum((q.').^(m.*(m+1)),2);
plot(m,cs);
grid on;

точность, которую дает нам сумма первых сто слагаемых на устраивает, но в дальнейших вычислениях нужно помнить, что этот способ не сработает для q в интервале от 0.99 до 1 (чем ближе к единице - тем больше слогаемых надо оставлять).

построим график от q

q = 0:.01:.99;
m = 0:100;
s = sum(q.^(m.'.*(m.'+1)),1);
plot(q,s);
grid on;

 

построим график для исходной функции со знаком и синусом, для разных значений i

q = 0:.1:.99;
m = 0:100;
N = 8;
for i = 0:N/2
    s = sum( ...
        ( (-1).^m.' .* q.^(m.'.*(m.'+1)) ) ...
        .* sin( (2*m.'+1)*pi*i/N )...
        ,1);
    %
    plot(q,s);
    hold on;
end
hold off;
grid on;
legend(compose('i = %d',0:N/2));

 

 

[C2000 2]

 

On 8/7/2021 at 10:55 AM, stealthisname said:

один из возможных способов расчета подобных выражений рассмотрим на примере суммы в числителе выражения

В принципе сумма в приведенном примере сходиться уже после 5 элементов, но результат почему-то некорректный.

 

 

Пробую рассчитать Elliptic Halfband Filter, по методике указанной в этой статье (страница 3, начиная с пункта (i)):

IEEE_Elliptic filter design for a class of generalized halfband filters.pdf

Беру за пример для расчёта Example II (стр.4) из этой же статьи, но коэффициенты никак не хотят совпадать с теми что автор рассчитал, и как следствие желаемые характеристики фильтра не получаются. С коэффициентами из статьи все гуд.

Прошу помощи у тех кто имеет больший опыт с матлабом и опыт в дизайне фильтров.

Возможно в статье ошибка? Например в формуле (36) по-моему должна быть сумма для m от 1 до inf, а не для i.

Вот мой скрипт, гляньте пожалуйста, что не так:

 

% (i) We are given input data
Ap = 0.001;     % is the maximum permissible passband magnitude deviation
As = 0.001;     % is the maximum permissible stopband magnitude deviation
Wp = 0.48*pi;   % frequency pass
Ws = pi - Wp;   % frequency stop
L = 6;          %

 

% (ii) Compute A as in (20).
A = min(As, (2*Ap - Ap^2)^0.5); % the required stopband magnitude deviation for the analog filter

 

%(iii) Compute k
k = tan(Wp/2)^2;
k_ = (1-k^2)^0.5;

 

%(iv) If
p = 0.5 * (1 - k_^0.5) / (1 + k_^0.5);
q = p + 2*(p^5) + 15*(p^9) + 150*(p^13);

% (iv) From (27) and (31), obtain N = 2L + 1;
N = 2*L + 1;

 

% (vi) Compute by approximating the expansion of (23):
ai = zeros(L,1);
for i = 1:L
syms m
f = ((-1)^m) * ( q^(m*(m+1)) ) * sin( ((2*m + 1)*pi*i) / N );
sum_1 = symsum(f, m, 0, 50);

syms m
f = ((-1)^m)* (q^(m^2)) * cos(2*m*pi*i/N);
sum_2 = symsum(f, m, 1, 50);

Om = (2*(q^0.25) * sum_1) / (1 + sum_2);

 

% (vii) Compute cos using (25) and (26) 
ri = abs( ( (1 - k * (Om^2))*(1 - (Om^2)/k) ) )^0.5;
cosO = ( ((-1)^(i+1)) * ri ) / (1 + Om^2);
ai_ = (1 - cosO) / (1 + cosO);

format long      %выводить много цифр после запятой
if (ai_ <= 1)
   ai(i) = ai_;
else
   ai(i) = 1/ai_;
end;

end

ai

 

[C2000 2]

Добрый день!

Хочется опять поднять тему. На all-pass фильтрах удалось получить нужную ФЧХ в диапазоне 40 - 4000Гц.

Но все же мысли о КИХ (фильтр Гильберта) не покидает.

Вот есть мысль, можно ли сделать перенос нужных частот подальше от 0, применить фильтр Гильберта а затем вернуть обратно? Тем самым можно было бы существенно снизить порядок этого самого фильтра. Да и для БИХ (на основе all-pass) это тоже актуально.

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

[petrov 6]

C2000

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

[C2000 2]

1 hour ago, petrov said:

C2000

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

Добрый день! Можно подробнее. Хотя есть сомнения что это возможно в реалтайме на МК (Cortex-M4 100MHz и это далеко не единственная задача которую он выполняет), с учётом что нужны все гармоники и интергармоники, и недопустим пропуск отсчётов.

Если даже мы получим DFT в реалтайме после каждого отсчёта с нужным шагом по частотам, то еще нужно для каждой выборки (после каждого отсчета) повернуть все гармоники на 90 градусов и рассчитать результат (суммарный) по всем частотам, а это совсем не проще с точки зрения нагрузки на ядро МК

[C2000 2]

А всё же что насчёт сдвига в частотной области поворота на 90 градусов и обратного восстановления? Это бесперспективно? Например сдвинуть на Fs/4 - по сути перенос "важных" частот в середину где порядок фильтра Гильберта будет на порядок:) меньше

[looser 8]

23 минуты назад, C2000 сказал:

А всё же что насчёт сдвига в частотной области поворота на 90 градусов и обратного восстановления? Это бесперспективно? Например сдвинуть на Fs/4 - по сути перенос "важных" частот в середину где порядок фильтра Гильберта будет на порядок:) меньше

Можно. Может вы все же внятно определите полосу частот и нужную точность измерения. Вообще-то гильберт 255-го порядка на м4 отъест примерно 4мгц производительности на 8кгц дискретизации.

[C2000 2]

32 minutes ago, looser said:

Можно. Может вы все же внятно определите полосу частот и нужную точность измерения. Вообще-то гильберт 255-го порядка на м4 отъест примерно 4мгц производительности на 8кгц дискретизации.

Важные частоты до 50 гармоники -> 2.5КГц, но с 2.5КГц то как раз проще, сложнее на 50Гц - 250Гц а они самые важные. Точность обусловленная ошибкой вычисления (есть естественно и другие) менее 0.1%. Вообще суммарная допустимая ошибка для реактивной мощности по ГОСТ 1% но хотелось бы чтобы в вычислениях она была минимальной, т.к. и без ошибки вычислений в большом динамическом диапазоне сигналов будет достаточная ошибка обусловленная аппаратными возможностями.

255 порядка явно не хватит чтобы не получить искажений по амплитуде на 50 - 250Гц 

32 minutes ago, looser said:

Можно. Может вы все же внятно определите полосу частот и нужную точность измерения. Вообще-то гильберт 255-го порядка на м4 отъест примерно 4мгц производительности на 8кгц дискретизации.

Как Вы это посчитали если не секрет? У вас получается что на каждое звено фильтра 2 такта МК, это не достижимо.  Нужно еще учитывать что 3 фазы, т.е. в любом случае нужно умножить на 3, а это лишь один из НЕМНОГИХ параметров который нужно рассчитать, не учитывая уже многих других задач помимо метрологии.

[petrov 6]

C2000

Можно подробнее.

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

https://electronix.ru/forum/index.php?app=forums&module=forums&controller=topic&id=23652&page=13#comment-1755048

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

_

Хотя есть сомнения что это возможно в реалтайме на МК (Cortex-M4 100MHz и это далеко не единственная задача которую он выполняет), с учётом что нужны все гармоники и интергармоники, и недопустим пропуск отсчётов.

Если даже мы получим DFT в реалтайме после каждого отсчёта с нужным шагом по частотам, то еще нужно для каждой выборки (после каждого отсчета) повернуть все гармоники на 90 градусов и рассчитать результат (суммарный) по всем частотам, а это совсем не проще с точки зрения нагрузки на ядро МК

У вас суперкомпьютер на котором можно сделать космический корабль и сесть на луну.