Сдвинуть сигнал на 90 градусов

[stealthisname 0]

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 0]

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

hilbert_Olli_Niemitalo (90 degree phse shift filter).pdf

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

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

 

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

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

[C2000 0]

 

....

[C2000 0]

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 0]

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 0]

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

 

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

[stealthisname 0]

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 0]

 

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