Jump to content

    
Sign in to follow this  
C2000

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

Recommended Posts

9 hours ago, C2000 said:

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

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

 

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

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

1880985709_.thumb.png.0549d8e52eba1f75b15b1de8303c9abd.png

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

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

 

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

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

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

k = -.9;

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

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

451162770_.thumb.png.ef4ba76568c2b0274ce4fc036d6a1c7c.png

 

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

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';

1468086986_.thumb.png.50b1ae71b9bade8a979c17042ae1ee70.png

 

Share this post


Link to post
Share on other sites

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

hilbert_Olli_Niemitalo (90 degree phse shift filter).pdf

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

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

 

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

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

Share this post


Link to post
Share on other sites
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 фильтром, но как, не понял пока что.

Share this post


Link to post
Share on other sites
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);

 

Share this post


Link to post
Share on other sites
54 minutes ago, C2000 said:

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

911319103_.thumb.jpg.0315ed7cb6ef0af05d7b94fc909a08e2.jpg

 

 

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

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

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

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

1562860502_.thumb.png.3c5e74c602ca7817a2cc75b3b79afc63.png

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

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

905852877_.thumb.png.51d337b75b6808b0ce7ed761d244c239.png

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

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

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

290370672_.thumb.png.da503758b2fe6abe459c161466672aa6.png

 

построим график для исходной функции со знаком и синусом, для разных значений 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));

1561984624_.thumb.png.2f85e7b0decace146224dd75bf0ecabf.png

 

 

Share this post


Link to post
Share on other sites

 

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

 

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this