Отправлено SM 10 июня 2004 г. 23:51

function tau=grp_del(b,a,w)
b=fliplr(b);
a=fliplr(a);
[b,a]=eqtflength(b,a);
Br=b;
Br(2:2:end)=0;
Br(3:4:end)=-Br(3:4:end);
Bi=b;
Bi(1:2:end)=0;
Bi(4:4:end)=-Br(4:4:end);
Ar=a;
Ar(2:2:end)=0;
Ar(3:4:end)=-Ar(3:4:end);
Ai=a;
Ai(1:2:end)=0;
Ai(4:4:end)=-Ar(4:4:end);
Ar=fliplr(Ar);
Ai=fliplr(Ai);
Br=fliplr(Br);
Bi=fliplr(Bi);
den_B=conv(Br,Br)+conv(Bi,Bi);
den_A=conv(Ar,Ar)+conv(Ai,Ai);
[num_B,tmp]=polyder(Bi,Br);
[num_A,tmp]=polyder(Ai,Ar);
tau=-polyval(num_B,w)./polyval(deb_B,w) + polyval(num_A,w)./polyval(deb_A,w);

• По любому - спасибо за алгоритм, было интересно — andy_P   (11.06.2004 17:09, пустое)
• Попробуйте grpdelay(b,a,N) - Она в Signal processing Toolbox входит — andy_P   (11.06.2004 16:29, пустое)
  • Не прокатит - он только для z-плоскости. Для аналоговых систем grpdelay в матлабе не работает. — SM   (11.06.2004 16:33, пустое)
    • Мне кажется, должно сработать - задержка по единичному кругу считается - поэтому что z-преобразование, что s, если у Вас АЧХ в дробно-рациональной форме есть... — andy_P   (11.06.2004 16:36, пустое)
      • Нэ работает! — SM   (11.06.2004 16:38, пустое)
        • Тьфу ты, соглаcен, там -pi->-Inf, +pi->+Inf. Преобразование требуется. Можно задать большое N и получить для любой разумной полосы, F<<fd. — andy_P   (11.06.2004 16:42, пустое)
        • Да и не может работать - в s-плоскости не по кольцу же считается.... — SM   (11.06.2004 16:41, пустое)
          • Точнее по кольцу, но оно не такой масштаб по частоте имеет — SM   (11.06.2004 16:48, пустое)
            • Там как раз билинейное преобразование работает - на чем я в DSP конфе наколося :-( - при f<<fd практически линейно, сильно нелинейно при больших f. Т.е. иммет смысл считать по небольшой части круга, где практически линейное преобразование. — andy_P   (11.06.2004 16:54, пустое)
              • Проще юзать специальные функции. Как то приведенная мной или freqs. Кстати (+) — SM   (11.06.2004 16:56, 115 байт)
                • С инвариантной ИХ в этом смысле еще проще - grpdelay нигде врать не будет (до fd пополам, разумеется), если при прототипировании алиасинга было мало(fd>>Fpass). На счет freqs() - согласен :-). Про прямые методы - как бы у нас система аналоговая :-) — andy_P   (11.06.2004 17:06, пустое)
  • Дело в том, что там алгоритм навороченный - он IIR сначала на секции второго порядка делит — andy_P   (11.06.2004 16:31, пустое)
• Небольшой глючок-с :) — SM   (11.06.2004 09:46, 84 байт)