stealthisname 7 6 августа, 2021 Опубликовано 6 августа, 2021 · Жалоба 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 3 6 августа, 2021 Опубликовано 6 августа, 2021 · Жалоба Нашел еще вот такую интересную статью: hilbert_Olli_Niemitalo (90 degree phse shift filter).pdf Вот только автор не рассказывает как он коэффициенты к этим фильтрам подбирал, только анализ результата. Есть ли какие то мысли, как эти фильтры рассчитывались? Там в конце дана формула одного звена: out(t) = a^2*(in(t) + out(t-2)) - in(t-2) Как получить АЧХ и ФЧХ именно по функции, в конечном итоге хотелось бы Си файл с функцией фильтра подцеплять и проверять что получается. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
C2000 3 6 августа, 2021 Опубликовано 6 августа, 2021 · Жалоба .... Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
C2000 3 6 августа, 2021 Опубликовано 6 августа, 2021 · Жалоба 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 6 августа, 2021 Опубликовано 6 августа, 2021 · Жалоба 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 3 7 августа, 2021 Опубликовано 7 августа, 2021 · Жалоба Подскажите как в Matlab посчитать вот это: При использовании symsum с аргументов inf Матлаб виснет. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
stealthisname 7 7 августа, 2021 Опубликовано 7 августа, 2021 · Жалоба 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 3 11 августа, 2021 Опубликовано 11 августа, 2021 · Жалоба 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 3 16 декабря, 2022 Опубликовано 16 декабря, 2022 · Жалоба Добрый день! Хочется опять поднять тему. На all-pass фильтрах удалось получить нужную ФЧХ в диапазоне 40 - 4000Гц. Но все же мысли о КИХ (фильтр Гильберта) не покидает. Вот есть мысль, можно ли сделать перенос нужных частот подальше от 0, применить фильтр Гильберта а затем вернуть обратно? Тем самым можно было бы существенно снизить порядок этого самого фильтра. Да и для БИХ (на основе all-pass) это тоже актуально. При переносе частот возникает проблема с зеркальным отображением относительно частоты на которую переносим. Да и сохраниться ли сдвиг 90 градусов будет после такого трюка, даже при наличии только одной (фундаментальной) частоты. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
petrov 7 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба C2000 Эта задача гораздо проще решается. Фильтры делаются в виде скользящего DFT, одно комплексное умножение на фильтр на входной отсчёт, вне зависимости от длины фильтров, фильтры вычисляются только в окрестности нужных гармоник, по выходам трёх фильтров основной гармоники считается её частота, ну и всё остально тривиально, любые соотношения гармоник, фазы комплексных сигналов и т. п. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
C2000 3 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба 1 hour ago, petrov said: C2000 Эта задача гораздо проще решается. Фильтры делаются в виде скользящего DFT, одно комплексное умножение на фильтр на входной отсчёт, вне зависимости от длины фильтров, фильтры вычисляются только в окрестности нужных гармоник, по выходам трёх фильтров основной гармоники считается её частота, ну и всё остально тривиально, любые соотношения гармоник, фазы комплексных сигналов и т. п. Добрый день! Можно подробнее. Хотя есть сомнения что это возможно в реалтайме на МК (Cortex-M4 100MHz и это далеко не единственная задача которую он выполняет), с учётом что нужны все гармоники и интергармоники, и недопустим пропуск отсчётов. Если даже мы получим DFT в реалтайме после каждого отсчёта с нужным шагом по частотам, то еще нужно для каждой выборки (после каждого отсчета) повернуть все гармоники на 90 градусов и рассчитать результат (суммарный) по всем частотам, а это совсем не проще с точки зрения нагрузки на ядро МК Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
C2000 3 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба А всё же что насчёт сдвига в частотной области поворота на 90 градусов и обратного восстановления? Это бесперспективно? Например сдвинуть на Fs/4 - по сути перенос "важных" частот в середину где порядок фильтра Гильберта будет на порядок:) меньше Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
looser 8 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба 23 минуты назад, C2000 сказал: А всё же что насчёт сдвига в частотной области поворота на 90 градусов и обратного восстановления? Это бесперспективно? Например сдвинуть на Fs/4 - по сути перенос "важных" частот в середину где порядок фильтра Гильберта будет на порядок:) меньше Можно. Может вы все же внятно определите полосу частот и нужную точность измерения. Вообще-то гильберт 255-го порядка на м4 отъест примерно 4мгц производительности на 8кгц дискретизации. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
C2000 3 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба 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 7 17 декабря, 2022 Опубликовано 17 декабря, 2022 · Жалоба 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 градусов и рассчитать результат (суммарный) по всем частотам, а это совсем не проще с точки зрения нагрузки на ядро МК У вас суперкомпьютер на котором можно сделать космический корабль и сесть на луну. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться