Jump to content

    

stealthisname

Участник
  • Content Count

    28
  • Joined

  • Last visited

Posts posted by stealthisname


  1. 10.10.2021 в 08:37, p_v сказал:

    Как определять "отсутствие" оборотов (или когда ниже допустимых, что актуально для старта).

    для того, чтобы определить отсутствие оборотов, можно сравнивать значение найденного максимума с порогом. порог для простоты можно рассчитывать по усредненным значениям fft.

    в качестве примера можно попробовать такой алгоритм на сигнале из Вашего файла hilda_15625Hz_zero_to_middle.

    я вычислил fft для 512 отсчетов на протяжении всего файла. для каждого результата убрал первые несколько отсчетов - в них мы ожидаем сильную низкочастотную помеху.

    получилась такая спектрограмма.

    1580207989_.thumb.png.f70c68946d30d54aa6c927843ea609e6.png356960761_.thumb.png.a9aa2b174ec1696571e3bdac07cd7b4a.png

    для каждого результата нашел среднее значение, и умножил на 12 - это порог, с которым сравнивается найденный максимум

    в начале файла шум не превышает порог

    1931207851_.thumb.png.0d0ae62d4ab920ba63d4ea52dc34dc0d.png

    во время включения уже виден нужный пик, но частота слишком быстро изменяется на интервале одного fft, поэтому порог не превышен

    609641312_.thumb.png.70bc36724e4e598312fb39bd00f1ff1f.png

    после выхода искомой частоты на постоянную все выборки дают максимум, превышающий порог

    1786396732_.thumb.png.98c8b8ab32419bbd9c095e5d62643910.png

    таким образом отсутствие оборотов не прошло алгоритм, PID/ADRC точно не подхватит ложный пик. Единственное, чем пришлось пожертвовать - скоростью захвата, так как мы определили включение немного позже действительного включения устройства.

     

    вот как выглядит спектрограмма с порогом

    977232755_.thumb.png.6aa997591994781a422f2a1692a02f22.png1495844636_.thumb.png.868aa84fada294cb76eb71285da40a31.png

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

    результат работы алгоритма - оценка частоты

    1368709200_.thumb.png.19440211771dcb8d2c55929d14b3334a.png198641435_.thumb.png.baf205553ff7b5346e5b7bfa9c3795cd.png

    аналогичные результаты для сигнала из файла hilda_15625Hz_zero_to_low  с тем же размером fft и с тем же множителем при вычислении порога

    143583384_.thumb.png.18c578ae2674ee1bab858125aeee4448.png2037838990_.thumb.png.a9748a6a60273be8257f44d8b7be9682.png

     

    228477947_.thumb.png.76fdfa9e7448364100e7db67659998ef.png163997626_.thumb.png.0efed0f06dc9ff6673f459ece8586b6f.png

  2. 10.10.2021 в 08:37, p_v сказал:

    Что все-таки смотреть после FFT, реальную часть или корень из суммы квадратов?

    если FFT выдает комплексные значения, то учитывать нужно и действительную и мнимую часть.

    по описанию алгоритма Вы ищете максимум из значений, поэтому вместо корня из суммы квадратов действительной и мнимой части, можно использовать просто сумму квадратов действительной и мнимой части. При этом положение максимума останется прежним, но можно будет сэкономить на операции взятия корня.

  3. 8 минут назад, RobFPGA сказал:

    Фактически  X это А, Yэто B,  перенос это PCIN, Out это либо PCOUT (если к идет переносу) либо P , например 

    спасибо

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

    схема получается такая

    2126749371_.png.3fc003bcc86023a7c1e38839c5557323.png.d68950d789f97f9d97648d916eafaf3b.png

    не считая обозначений входов, все было правильно

  4. 19 минут назад, RobFPGA сказал:

    Ну так DSP может суммировать  3 числа: P = AB+C+PCIN каждое по 48 бит (AB тут не произведение а просто конкатенация входов).  

    на первой схеме у автора темы таких входов нет, есть X1 Y2 и перенос. Есть ли возможность суммировать числа в DSP так, как указано автором темы на первой схеме?

     

  5. Только что, RobFPGA сказал:

    Что значит "не укладываются входы DSP на последней стадии сложения"? 

    "Результат всей суммы укладывается в 36 бит - это по мат. модели. "

    на последней DSP при таком суммировании будет приходить два сигнал с разрядностью 35 бит, а в блоке DSP только один сигнал переноса, который может вместить такую разрядность.

    4 минуты назад, RobFPGA сказал:

    Удачи! Rob.

    Вам тоже, удачи!

  6. 1 час назад, alexPec сказал:

    Результат всей суммы укладывается в 36 бит - это по мат. модели.

    в схеме, которую Вы предлагаете реализовать, не укладываются входы DSP на последней стадии сложения

    2126749371_.png.3fc003bcc86023a7c1e38839c5557323.png.2cf489525b1afe8fc95a03dca9e9e600.png

    если перенос с достаточной разрядностью для суммирования в DSP всего один, то такая схема не реализуемая

     

  7. по условию задачи и по обсуждению выше, в цифровой обработке допускается использовать БИХ фильтры второго порядка (только если нет высоких требований к точности коэффициентов и к точности вычислений)

    я отфильтровал сигналы из hilda_50kHz_rpm_high и hilda_50kHz_rpm_low в ВЧ фильтре с коэффициентами 

    b = [1,-1]; a = [1,-0.9375]; 

    с такой АЧХ:

    23939753_.thumb.png.a7dbf642e016e9510092a184c3e60607.png

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

    Скрытый текст
    
    f =   400; b = [  0.03281759059,0, -0.03281759059]; a = [1,     -1.9340458,         0.9375];
    f =   900; b = [  0.03156366737,0, -0.03156366737]; a = [1,   -1.924120029,         0.9375];
    f =  1400; b = [   0.0313784307,0,  -0.0313784307]; a = [1,   -1.906600635,         0.9375];
    f =  1900; b = [  0.03131856954,0, -0.03131856954]; a = [1,    -1.88155676,         0.9375];
    f =  2400; b = [  0.03129199717,0, -0.03129199717]; a = [1,    -1.84908724,         0.9375];
    f =  2900; b = [  0.03127793598,0, -0.03127793598]; a = [1,   -1.809320217,         0.9375];
    f =  3400; b = [  0.03126961201,0, -0.03126961201]; a = [1,   -1.762412634,         0.9375];
    f =  3900; b = [  0.03126428538,0, -0.03126428538]; a = [1,   -1.708549613,         0.9375];
    f =  4400; b = [  0.03126067664,0, -0.03126067664]; a = [1,   -1.647943728,         0.9375];
    f =  4900; b = [  0.03125812333,0, -0.03125812333]; a = [1,   -1.580834161,         0.9375];

     

    АЧХ фильтров:

    2051006057_.thumb.png.a6e69bee6058ef93017062b1fb9262e5.png

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

    процесс такого накопления в hilda_50kHz_rpm_high:

    1770714601_.thumb.png.1f473eecfa47ae2c3ac88fdcca63d375.png

    результат такого накопления:

    17661489_.thumb.png.29edf940a7a0bf9a1afc1ff101023547.png

    результат накопления для hilda_50kHz_rpm_low:

    1454655599_.thumb.png.c6308a02e1db58fef6efe5c7704a7507.png

     

    в сигнале hilda_50kHz_rpm_high больше всего накопил фильтр на частоте 3900. в сигнале hilda_50kHz_rpm_low больше всего накопил фильтр на частоте 400.

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

    если потребуется больше точность вычисления и будет позволять скорость вычислений, можно увеличить количество фильтров. больше фильтров - больше точность вычисления.

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

     

  8. 14 часов назад, p_v сказал:

    Кому интересно, исходный сигнал - ток коллекторного мотора, по которому пытаемся определить обороты (в репозитории есть дампы). При вращении происходит перекоммутация обмоток, и это вызывает весьма заметные колебания, пропорциональные оборотам и количеству полюсов. Из помех - какие-то адские рандомные "прострелы", плюс 100 герц из сети.

     

    Здравствуйте!

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

    я скачал файл hilda_50kHz_rpm_high.txt и попробывал применить медианный фильтр.

    оригинальный сигнал:

    1154958688_.thumb.png.9442dd3f01df5fe04cca106ccbf9656f.png

    результат обработки медианным фильтром третьего порядка:

    1960551823_.thumb.png.c88b869528b03390c46fd32cf103427c.png

    результат обработки медианным фильтром пятого порядка:

    1517149767_.thumb.png.cfff8bbbbbf79c567a81e380461f27f5.png

    по результатам обработки видно, что тут можно рассмотреть применение медианного фильтра до фильтрации ФНЧ и ФВЧ.

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

     

  9. 18 минут назад, turnip сказал:

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

    если я правильно понимаю, по описанию работы шины генератор импульсов, счетчик и дешифратор (левый верхний в схеме) не могут переключаться во время передачи, так как иначе следующее в цепочке переключений устройство будет считать, что оно может передавать данные по  общей шине. Как подобные коллизии решаются, если обратной связи устройств с генератором импульсов по схеме нет? что помешает счетчику переключиться во время передачи данных?

  10. 31 минуту назад, Александр77 сказал:

    этот вывод громоздок и попытка рассмотреть небольшую область проблематична.

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

    s = 0:1e-2:1;
    t = 0:1e-2:100;
    z = sin(pi*s).*cos(pi*t.'+2*pi*10*s)+cos(pi*t.'/10)/2;
    surf(s,t,z,'EdgeAlpha',0);
    view(0,90);

    735240660_.thumb.png.07546ea7bf5f4c20f9da17ea294d2b13.png

    разглядеть детально небольшую область в таком масштабе действительно проблематично, поэтому выбираем Zoom In, Vertical Zoom

    1093089070_.thumb.png.66bf82558dcf2f7b4d476ee02134eb14.png

    любой участок спектрограммы приближается в пару кликов, все изменения графика вдоль оси времени - отлично отслеживаются

    563372970_.thumb.png.46da2a8f8047f7e20d2803ccd2c6a571.png

    для плавного скольжения по спектрограмме вдоль оси времени выбираем Pan, Vertical Pan

    167140888_.thumb.png.4e251998b77dab157f4a22e223a690b3.png

    довольно удобный способ исследования спектрограмм

  11. 32 минуты назад, Александр77 сказал:

    В матлабе есть трехмерные графики интенсивности (mesh, surf), с цветовой палитрой. При больших объемах данных, эти графики неповоротливы и плохо масштабируемы (на мой взгляд).

    Меня интересует есть ли возможность отображения интенсивности в двумерном формате?

    если уже есть вычисление данных и построение с помощью surf, например такое

    t = 0:1e-2:1;
    z = sin(2*pi*t).*cos(pi*t.');
    surf(z);

    937718821_.thumb.png.c90e547646d9f47ad4ebb52d413079b9.png

    то получить двухмерный вариант можно например командой view(0,90), вот так

    t = 0:1e-2:1;
    z = sin(2*pi*t).*cos(pi*t.');
    surf(z);
    view(0,90);

    1831157109_.thumb.png.7a67f810dd3f844dcef5cc8bdb925f57.png

    таким образом будет выведен двухмерный вариант графика интенсивности, но так же сохранится возможность посмотреть и трехмерный график, с помощью кнопки Rotate 3D

     

     

  12. 2 hours ago, _sda said:

    А попробуйте сгенерировать HDL для коэффициента не 7/5 а 13/5. У меня получилось вот так:

    для проверки преобразования частоты с коэффициентом 7/5 я использовал импульсную характеристику по такому образцу

    b  = firls(N, [0, 0.9*(1/7), 1.1*(1/7), 1], [1 1 0 0], [Wpass Wstop]);

    специально, чтобы фильтр оставлял примерно 1/7 полосы после интерполирования в 7 раз

    для преобразования частоты с коэффициентом 13/5 такая ИХ не подойдет, так как необходимо , чтобы фильтр оставлял примерно 1/13, с ИХ по такому образцу

    b  = firls(N, [0, 0.9*(1/13), 1.1*(1/13), 1], [1 1 0 0], [Wpass Wstop]);

    при этом для достижения заданных характеристик фильтра потребуется больший порядок (N)

  13. удалось сгенерировать  похожий в матлабе в Filter Designer

    после экспорта в .m файл генерируется с такого скрипта

     

    % MATLAB Code
    % Generated by MATLAB(R) 9.2 and the DSP System Toolbox 9.4.
    % Generated on: 25-Aug-2021 23:28:00
    
    % FIR least-squares Lowpass filter designed using the FIRLS function.
    
    % All frequency values are in Hz.
    Fs = 14;  % Sampling Frequency
    
    N     = 70;   % Order
    Fpass = 0.9;  % Passband Frequency
    Fstop = 1.1;  % Stopband Frequency
    Wpass = 1;    % Passband Weight
    Wstop = 1;    % Stopband Weight
    
    % Calculate the coefficients using the FIRLS function.
    b  = firls(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop]);
    intf = 7;                     % Interpolation Factor
    decf = 5;                     % Decimation Factor
    
    Hd = dsp.FIRRateConverter( ...
        'InterpolationFactor', intf, ...
        'DecimationFactor', decf, ...
        'Numerator', b, ...
        'CoefficientsDataType', 'Custom', ...
        'CustomCoefficientsDataType', numerictype([],16,17));
    
    generatehdl(Hd, 'ResetType', 'Synchronous',... 
                   'EDAScriptGeneration', 'off',... 
                   'Name', 'Filter_fract_rate_fdatool',... 
                   'TargetLanguage', 'Verilog',...
                   'InputDataType',numerictype(1,16,15));
    

     

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

    454622929_.thumb.png.2eb84ae11486ed949ff278735d52272d.png

    2028842182_.thumb.png.a250e6144a4bf699682b1e50bf5a9f62.png

    в самом Filter Designer были выбраны такие настройки

    1740877225_.thumb.png.f665b61b15afb2e6789c94b32f4627fe.png

    1704000282_.thumb.png.0c3b747d69961a72d5d43a8ce725e6b8.png

    1740883322_.thumb.png.830269be0ae89236605cb5c651267f55.png

     

     

     

  14. 6 hours ago, _sda said:

    Кто нибудь пользовался сгенерёнными в матлабе HDL-файлами? Как впечатление?

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

     

    очень похоже, что для Вашего случая подойдет блок FIR Rate Conversion HDL Optimized

    для проверки блока с дискретизацией 7/5 я использовал такие параметры

    706177632_.png.2afbfcf5e2c51ac54d5f262ef7239be1.png

     

    результаты моделирования с входной дискретизацией как в Вашем случае (входной отсчет раз в два такта)

    1897780915_.thumb.png.9ca6d3c3073911c592ff886cc6c1d1d3.png

    930697728_.thumb.png.d1c1025215720d825bf27e3ac1576fcc.png

     

     

     

  15. 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

     

     

  16. 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);
    
    

     

  17. 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

     

  18. 6 hours ago, C2000 said:

    Можно ли также построить АЧХ и ФЧХ двух фильтров каскадом или разницу между выходами двух фильтров.

    ИХ двух фильтров каскадом - свертка ИХ отдельных фильтров. по рассчитанной ИХ можно построить АЧХ и ФЧХ

    для КИХ фильтров в матлабе можно рассчитать например так

    % b1, b2 - ИХ отдельных фильтров в каскаде
    
    % общая ИХ
    b = conv(b1,b2);
    
    % вывод
    fvtool(b1,1,b2,1,b,1);

     

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

    % b1, b2 - коэфф числителя передаточной
    % a1, a2 - коэфф знаменателя передаточной
    
    % общая передаточная
    b = conv(b1,b2);
    a = conv(a1,a2);
    
    % вывод
    fvtool(b1,a1,b2,a2,b,a);

     

    ИХ суммы/разности сигналов с двух фильтров - это сумма/разность ИХ отдельных фильтров. По рассчитанной ИХ можно построить АЧХ и ФЧХ

  19. 6 hours ago, C2000 said:

    Можно код для Matlab в качестве примера, где сначала рассчитывается фильтр, потом к нему добавляется задержка в сэмплах и это всё выводиться.

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

    % FIR least-squares Hilbert transformer
    % filter designed using the FIRLS function.
    
    Fs = 8;  % Частота дискретизации (кГц)
    
    N =  40; % порядок фильтра
    F = [0.005, Fs/2-0.005]; % Frequency Vector
    A = [1, 1];              % Amplitude Vector
    W = [1];                 % Weight Vector
    
    b  = firls(N, F/(Fs/2), A, W, 'hilbert');
    
    %
    b0 = zeros(size(b));
    b0((end+1)/2) = 1;
    hfvt = fvtool(b0,1,b,1);
    hfvt.Analysis = 'freq';
    hfvt.PhaseUnits = 'degrees';
    
    
    

    тут b - импульсная характеристика рассчитанного фильтра, b0 - ИХ задержки на половину порядка рассчитанного фильтра.

    половина коэффициентов действительно получается нулевыми

    1353595411_.thumb.png.24a760ec9d84998e44914be7e0ac0c36.png

     

    я не нашел, как в fvtool построить отношение амплитудных характеристик двух фильтров и разность фазовых характеристик

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

    w = 0:1e-4:.5; % цифровая нормированная частота
    f = w*Fs; % частота (кГц)
    
    % комплексные частотные харатеристики
    h  = sum(b .'.*exp(-1i*2*pi*w.*(1:N+1).'));
    h0 = sum(b0.'.*exp(-1i*2*pi*w.*(1:N+1).'));
    
    % АЧХ (дБ)
    h_magn = mag2db(abs(h));
    h0_magn = mag2db(abs(h0));
    
    % ФЧХ (°)
    h_phas = rad2deg(unwrap(angle(h)));
    h0_phas = rad2deg(unwrap(angle(h0)));
    
    %
    % plot(f,h0_magn,f,h_magn);
    % plot(f,h0_phas,f,h_phas);
    
    % разность ФЧХ
    plot(f,h0_phas-h_phas);
    

     

  20. 11 hours ago, C2000 said:

    Нужно посчитать мощности в сети.

    Активная просто, интеграл от U*I

    А вот для расчёта реактивной напряжение нужно сдвинуть на 90 град, причем в идеале нужно подвинуть и гармоники.

     

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

     

    построим модель сигнала. Сначала сгенерим на нужной тактовой частоте нужную основную гармонику. Смоделируем вторую с заданной разностью фаз. (здесь я беру медленно линейно возрастающую до 60°).

    эти два сигнала будут имитировать принятые напряжение и ток

    Spoiler

    res_00.thumb.png.7ae4429addb59a6bcb8069f4ab91e87b.png

     

    добавим дополнительных кратных гармоник, добавим шум.

     

    Spoiler

    res_01.thumb.png.e4539bc0a1f74363c2b2bb2cbf6497ad.png

    res_02.png.56c28329360ba5febaf777a067a93694.png

     

    в обработке сигнала перенесем с заданной известной частоты в ноль, при этом получится комплексный сигнал

    Spoiler

    res_03.thumb.png.9fe80b2adfdf1433fe60b70feba745e7.png

    res_04.png.b4a201b2e9f9f655ff6be22f71ef7df3.png

     

    результат продецимируем, например в CIC фильтре

    Spoiler

    res_05.thumb.png.f7a005be724575c896f20a9951d09f94.png

    res_06.png.37539a999612f4ddd7fce8b6a865287b.png

     

    по отфильтрованному продецимированному сигналу, найдем произведение напряжения и комплексно сопряженного с током

    Spoiler

    res_07.thumb.png.725c23c09707aafa6614f76b004f21c7.png

    res_08.png.2deb94125bb516a125f883eb8237b993.png

     

    накопим значения для очищения от шума (можно с помощью еще одно CIC фильтра)

    Spoiler

    res_09.thumb.png.08fb08b42ee30bf0366aef50ce99f60c.png

    res_10.png.e1a216b754d73a34dd480f7c610c5c4e.png

     

    эта комплексная величина как-раз пропорциональна мощности гармоники: действительная часть пропорциональна активной составляющей, мнимая часть пропорциональна реактивной составляющей.

    теперь можем рассчитать значение коэффициента, на который нужно домножить активную мощность, чтобы получить реактивную

    делим

    Spoiler

    res_11.thumb.png.25cfbb392d0f579cf6e0636e8e1ee769.png

    res_12.png.e2f2b6fa3be0b47c787643223489d0a2.png

     

    чтобы убедиться, что алгоритм расчета устойчив к изменению входной частоты, изменим входную частоту на 58 Гц

     

    Spoiler

    res_13.thumb.png.e7192ae722dff9d234e74502c04680f1.png

    res_14.thumb.png.3d186df216f83b819620d044093f0b14.png

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

     

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

     

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

  21. 10 hours ago, C2000 said:

    Смотрел на фильтр гильберта в Матлабе, но так и не понял как его прикрутить, и вообще то ли это.

    фильтр гильберта в Матлабе подходит для указанной задачи, но с некоторыми ограничениями.

    в задаче указано, что максимум порядок фильтра - несколько десятков, поэтому расчитаем фильтр гильберта для предельного случая, когда порядок фильтра равен 90, больше мы точно брать не можем.

    частота дискретизации 8кГц. Полоса фильтра, с учетом порядка фильтра, получилась от 50 Гц до 2кГц.

    res_00.thumb.png.6bc078f91a2a4a4e211615313cc291db.png

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

    посмотрим ИХ

    Spoiler
    
    % FIR least-squares Hilbert transformer
    % filter designed using the FIRLS function.
    
    % All frequency values are in kHz.
    Fs = 8;  % Sampling Frequency
    
    N =  90;                            % Order
    F = [0, 0.0001, 0.005, 2, 2.1, 4];  % Frequency Vector
    A = [0, 0, 1, 1, 0, 0];             % Amplitude Vector
    W = [1, 10, 1];                     % Weight Vector
    
    b  = firls(N, F/(Fs/2), A, W, 'hilbert');
    
    %
    b0 = zeros(size(b));
    b0((end+1)/2) = 1;
    hfvt = fvtool(b0,1,b,1);
    hfvt.Analysis = 'impulse';

     

    res_01.thumb.png.1057e9514ab3d1bf8bc97a2ae4640e79.png

     

    построим АЧХ

    hfvt = fvtool(b0,1,b,1);
    hfvt.Fs = 8e3;
    hfvt.Analysis = 'magnitude';

    res_03.thumb.png.9c293f681f8f3dbbf22240cc4f3d54b9.png

    построим ФЧХ

    hfvt = fvtool(b0,1,b,1);
    hfvt.Fs = 8e3;
    hfvt.Analysis = 'phase';
    hfvt.PhaseUnits = 'degrees';

    res_04.thumb.png.859fdf4af4c465f6c214a6ae6b7a9287.png

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

     

    порядок фильтра равный 90 был взят как предельный возможный, при уменьшении порядка характеристики фильтра будут ухудшаться: сдвиг фазы все еще будет 90°, но полоса будет уменьшаться.

    для изучения влияния порядка фильтра на характеристики, сделаем аналогичные расчеты для более реального порядка фильтра, равного 30

     

    Spoiler
    
    % FIR least-squares Hilbert transformer
    % filter designed using the FIRLS function.
    
    % All frequency values are in kHz.
    Fs = 8;  % Sampling Frequency
    
    N =  30;                            % Order
    F = [0, 0.0001, 0.005, 2, 2.1, 4];  % Frequency Vector
    A = [0, 0, 1, 1, 0, 0];             % Amplitude Vector
    W = [1, 10, 1];                     % Weight Vector
    
    b  = firls(N, F/(Fs/2), A, W, 'hilbert');
    
    hfvt = fvtool(b0,1,b,1);
    hfvt.Fs = 8e3;
    hfvt.Analysis = 'magnitude';

    res_05.thumb.png.1f2df00c9a2aa661cb6ec55e43c53509.png

    
    hfvt = fvtool(b0,1,b,1);
    hfvt.Fs = 8e3;
    hfvt.Analysis = 'phase';
    hfvt.PhaseUnits = 'degrees';

    res_06.thumb.png.394c13b16043d0ba4a7191890407aa01.png

     

  22. On 6/7/2021 at 6:21 AM, Nik_Su said:

    Все сделал, даже получилось выделить огибающую, а вот с восстановлением сигнала ни как не получается, пытался применить формулы atan2(Q,I), diff(atan(Q./I)), diff(atan2(Q, I)) ничего не получается.

    для вычисленных Вами квадратур формула atan2(Q,I) применяется и дает возможность увидеть исходную фазовую модуляцию

    я применял формулу к такому участку вычисленных Вами квадратур

    plot(t,I_signal,t,Q_signal);

    res_00.thumb.png.a61df8430366af758f0c2fb6792485f0.png

    закон изменения фазы принятого сигнала

    plot(atan2(I_signal,Q_signal));

     

    res_01.thumb.png.c3fd3676bb4ad5420372844a937aab9c.png

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

    plot(unwrap(atan2(I_signal,Q_signal)));

     

    res_02.thumb.png.7a4b3eaa636a3bf5f99d3846df2f0baf.png

    закон изменения фазы - порабола, как и положено для ЛЧМ сигнала. Также убедимся, что закон изменения частоты линейный

    plot(diff(unwrap(atan2(I_signal,Q_signal))));

    res_03.thumb.png.2323ce9c7719e787144175390aa4d59b.png

    по результатам расчетов видно, как именно использованный ранее в расчетах НЧ фильтр искажает закон изменения частоты сигнала при таких вычислениях

  23. 12 hours ago, Halfback said:

    Пока мне не очень понятно насколько корректно будет работать описание внутри posegde и negedge если их засунуть в один posegde. Надо проверять, т.к. вычисление должно происходить за один такт.

    надо проверять, и в данном случае можно разобраться на примере вашего фильтра

    тестбенч похож на предыдущий, но упрощен входной сигнал

    Spoiler
    
    `timescale 1 ns / 1 ps
    
    module MovAverageFilter_tb();
    
    reg clk = 0;
    
    always #(10.0/2)
    begin
        clk = ~clk;
    end
    
    reg [15:0] datain;
    
    
    MovAverageFilter 
    #(
        .FIFO_LEN_LOG(4)//,
    )
    MovAverageFilter_i
    (
    	.clk(clk),
    	.datain(datain)//,
    );
    
    MovAverageFilter_ram 
    #(
        .FIFO_LEN_LOG(4)//,
    )
    MovAverageFilter_ram_i
    (
    	.clk(clk),
    	.datain(datain)//,
    );
    
    initial
    begin
        datain = 0;
    end
    
    always @(posedge clk)
    begin
        
        datain = 1;
        
    end
    
    endmodule

     

    сначала посмотрим вариант с такой конструкцией

        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[k];        
        end
        
        always @(negedge clk)
        begin
            FIFO[k] <= datain;
            k <= k + 1;
        end
    Spoiler
    
     module MovAverageFilter_ram
    #(
        parameter FIFO_LEN_LOG = 4//,
    )
    (
        input  clk,
        input  [15:0] datain,
        output [15:0] dataout//,
    );
        localparam FIFO_LEN = (1 << FIFO_LEN_LOG);
        
        reg [15:0] FIFO [0:FIFO_LEN-1];
        reg [FIFO_LEN_LOG+15:0] fifo_sum = 0;
        integer i;
         
        reg [FIFO_LEN_LOG-1:0] k = 0;
        
        initial
        begin
            for(i = FIFO_LEN-1; i >= 0; i=i-1)  FIFO[i] = 16'd0;
        end
        
        
        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[k];        
        end
        
        always @(negedge clk)
        begin
            FIFO[k] <= datain;
            k <= k + 1;
        end
        
        wire [15:0] FIFO_k;
        assign FIFO_k = FIFO[k];
        
        assign dataout = fifo_sum >> FIFO_LEN_LOG;
        
    endmodule 

     

    фильтр работает не так как задумывалось

    res_01.thumb.png.97967acaf237774412c52880069781c4.png

    накопление в фильтре идет, но размер окна равен 15, а не 16. счетчик (k) расчитывается в @(negedge clk), а используется как в @(negedge clk), так и в @(posedge clk). задержка в памяти с учетом этого получается на такт меньше, и тут для получения нужной задержки требовалось бы еще задержать счетчик (k) до  @(posedge clk) для использования в чтении с памяти, или брать предыдущий отсчет считывая (k-1) или задерживать выход памяти, или делать что-то еще, но

    но посмотрим, что будет с конструкцией

        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[k];        
        end
        
        always @(posedge clk)
        begin
            FIFO[k] <= datain;
            k <= k + 1;
        end
    
    Spoiler
    
     module MovAverageFilter_ram
    #(
        parameter FIFO_LEN_LOG = 4//,
    )
    (
        input  clk,
        input  [15:0] datain,
        output [15:0] dataout//,
    );
        localparam FIFO_LEN = (1 << FIFO_LEN_LOG);
        
        reg [15:0] FIFO [0:FIFO_LEN-1];
        reg [FIFO_LEN_LOG+15:0] fifo_sum = 0;
        integer i;
         
        reg [FIFO_LEN_LOG-1:0] k = 0;
        
        initial
        begin
            for(i = FIFO_LEN-1; i >= 0; i=i-1)  FIFO[i] = 16'd0;
        end
        
        
        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[k];        
        end
        
        always @(posedge clk)
        begin
            FIFO[k] <= datain;
            k <= k + 1;
        end
        
        wire [15:0] FIFO_k;
        assign FIFO_k = FIFO[k];
        
        assign dataout = fifo_sum >> FIFO_LEN_LOG;
        
    endmodule 

     

    все работает так как предполагалось изначально

    res_02.thumb.png.1e08d1c2f4f07df2b70056fb2b67f445.png

    счетчик (k) расчитывается  в @(posedge clk) и используется в @(posedge clk), кольцевой буфер дает задержку ровно на 16 тактов, как и положено, размер окна равен 16

  24. по заданию требуется возможность генерации модуля с разным размером окна и разрешены размеры окна по степеням двойки, поэтому в модуле можно в качестве параметра задавать степень (FIFO_LEN_LOG) и на основе ее рассчитывать размер окна (FIFO_LEN)

        parameter FIFO_LEN_LOG = 4;
        localparam FIFO_LEN = (1 << FIFO_LEN_LOG);
    

    обьявление, инициализация и работа внутренней памяти (FIFO) останется прежней, но с учетом новых параметров

        reg [15:0] FIFO [0:FIFO_LEN-1];
        integer i;
    
        initial
        begin
            for(i = FIFO_LEN-1; i >= 0; i=i-1)
            begin
                FIFO[i] = 16'd0;
            end
        end
        
        always @(posedge clk)
        begin
    
            for(i = FIFO_LEN-1; i > 0; i=i-1)
            begin
                FIFO[i] <= FIFO[i-1];
            end
            
            FIFO[0] <= datain;
            
        end

    как было замечено выше в обсуждении темы, для вычисления результата суммирования (fifo_sum) с использованием предыдущего результата рассчета количество сложений может быть уменьшено

        reg [FIFO_LEN_LOG+15:0] fifo_sum = 0;
        
        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[FIFO_LEN-1];    
        end

    деление результата суммирования также можно записать с учетом параметров

        assign dataout = fifo_sum >> FIFO_LEN_LOG;

    весь модуль у меня получился таким

    Spoiler


    
    module MovAverageFilter
    #(
        parameter FIFO_LEN_LOG = 4//,
    )
    (
        input  clk,
        input  [15:0] datain,
        output [15:0] dataout//,
    );
        localparam FIFO_LEN = (1 << FIFO_LEN_LOG);
        
        reg [15:0] FIFO [0:FIFO_LEN-1];
        integer i;
    
        initial
        begin
            for(i = FIFO_LEN-1; i >= 0; i=i-1)
            begin
                FIFO[i] = 16'd0;
            end
        end
        
        always @(posedge clk)
        begin
    
            for(i = FIFO_LEN-1; i > 0; i=i-1)
            begin
                FIFO[i] <= FIFO[i-1];
            end
            
            FIFO[0] <= datain;
            
        end
        
        reg [FIFO_LEN_LOG+15:0] fifo_sum = 0;
        
        always @(posedge clk)
        begin
            fifo_sum <= fifo_sum + datain - FIFO[FIFO_LEN-1];	
        end
        
        assign dataout = fifo_sum >> FIFO_LEN_LOG;
        
    endmodule
    


     

    модуль проверялся вот таким тестбенчем

    Spoiler
    
    `timescale 1 ns / 1 ps
    
    module MovAverageFilter_tb();
    
    reg clk = 0;
    
    always #(10.0/2)
    begin
        clk = ~clk;
    end
    
    reg [15:0] datain;
    
    
    MovAverageFilter_16 MovAverageFilter_16_i
    (
    	.clk(clk),
    	.datain(datain)//,
    );
    
    MovAverageFilter 
    #(
        .FIFO_LEN_LOG(4)//,
    )
    MovAverageFilter_i_4
    (
    	.clk(clk),
    	.datain(datain)//,
    );
    
    MovAverageFilter 
    #(
        .FIFO_LEN_LOG(6)//,
    )
    MovAverageFilter_i_6
    (
    	.clk(clk),
    	.datain(datain)//,
    );
    
    real signal;
    real signal_noise;
    
    
    initial
    begin
        datain = 0;
        signal = 0;
        signal_noise = 0;
    end
    
    always @(posedge clk)
    begin
        
        //
        signal = (1-1e-5)*signal + (5e-8)*$random();
        signal_noise = signal + (1e-6)*$random() + 2.0**15;
        signal_noise = signal_noise + (1e-6)*$random();
        
        //
        datain = signal_noise;
        
        if (signal_noise>=(2.0**16-1))
        begin
            datain = (2.0**16-1);
        end
        if (signal_noise<=(0))
        begin
            datain = (0);
        end
        
    end
    
    endmodule

    тут

    (signal) - низкочастотный сигнал
    (signal_noise) - зашумленный низкочастотный сигнал

    (MovAverageFilter_16_i.dataout) - выход оригинального блока
    (MovAverageFilter_i_4.dataout) - выход параметризованного варианта с размером окна как у оригинального
    (MovAverageFilter_i_6.dataout) - выход параметризованного варианта с большим размером окна

    (MovAverageFilter_16) - оригинальный блок для сравнения результатов

    Spoiler
    
    
    module MovAverageFilter_16
    (
        input  clk,
        input  [15:0] datain,
        output [15:0] dataout//,
    );
        
        reg [15:0] FIFO [0:15];
        integer i;
    
        initial
        begin
            for(i = 15; i >= 0; i=i-1)
            begin
                FIFO[i] = 16'd0;
            end
        end
        
        always @(posedge clk)
        begin
    
            for(i = 15; i > 0; i=i-1)
            begin
                FIFO[i] <= FIFO[i-1];
            end
            
            FIFO[0] <= datain;
            
        end
        
        reg [19:0] fifo_sum = 0;
        
        always @(posedge clk)
        begin
            fifo_sum <=
            (
                +FIFO[ 0]+FIFO[ 1]+FIFO[ 2]+FIFO[ 3]
                +FIFO[ 4]+FIFO[ 5]+FIFO[ 6]+FIFO[ 7]
                +FIFO[ 8]+FIFO[ 9]+FIFO[10]+FIFO[11]
                +FIFO[12]+FIFO[13]+FIFO[14]+FIFO[15]
            );	
        end
        
        assign dataout = fifo_sum >> 4;
        
    endmodule

     

     

    результаты моделирования

    res_00.thumb.png.17e00d3146d97c30fb58654536ff16d6.png

     

    сигнал с выхода параметризованного варианта с размером окна как у оригинального повторяет сигнал с выхода оригинального блока, сигнал с выхода параметризованного варианта с большим размером окна отсеивает больше шума