Jump to content

    

stealthisname

Участник
  • Content Count

    28
  • Joined

  • Last visited

Everything posted by stealthisname


  1. для того, чтобы определить отсутствие оборотов, можно сравнивать значение найденного максимума с порогом. порог для простоты можно рассчитывать по усредненным значениям fft. в качестве примера можно попробовать такой алгоритм на сигнале из Вашего файла hilda_15625Hz_zero_to_middle. я вычислил fft для 512 отсчетов на протяжении всего файла. для каждого результата убрал первые несколько отсчетов - в них мы ожидаем сильную низкочастотную помеху. получилась такая спектрограмма. для каждого результата нашел среднее значение, и умножил на 12 - это порог, с которым сравнивается найденный максимум в начале файла шум не превышает порог во время включения уже виден нужный пик, но частота слишком быстро изменяется на интервале одного fft, поэтому порог не превышен после выхода искомой частоты на постоянную все выборки дают максимум, превышающий порог таким образом отсутствие оборотов не прошло алгоритм, PID/ADRC точно не подхватит ложный пик. Единственное, чем пришлось пожертвовать - скоростью захвата, так как мы определили включение немного позже действительного включения устройства. вот как выглядит спектрограмма с порогом видно, что алгоритм четко определяет только нужные пики искомой частоты и срабатывает только при включенном устройстве. результат работы алгоритма - оценка частоты аналогичные результаты для сигнала из файла hilda_15625Hz_zero_to_low с тем же размером fft и с тем же множителем при вычислении порога
  2. если FFT выдает комплексные значения, то учитывать нужно и действительную и мнимую часть. по описанию алгоритма Вы ищете максимум из значений, поэтому вместо корня из суммы квадратов действительной и мнимой части, можно использовать просто сумму квадратов действительной и мнимой части. При этом положение максимума останется прежним, но можно будет сэкономить на операции взятия корня.
  3. спасибо я ошибся, сказав, что такая сумма не возможна на DSP. вычисления у автора темы могут быть реализованы так как указано на первой схеме. но в этой схеме есть незначительная неточность, которая может из-за невнимательности ввести в заблуждение. схема получается такая не считая обозначений входов, все было правильно
  4. на первой схеме у автора темы таких входов нет, есть X1 Y2 и перенос. Есть ли возможность суммировать числа в DSP так, как указано автором темы на первой схеме?
  5. "Результат всей суммы укладывается в 36 бит - это по мат. модели. " на последней DSP при таком суммировании будет приходить два сигнал с разрядностью 35 бит, а в блоке DSP только один сигнал переноса, который может вместить такую разрядность. Вам тоже, удачи!
  6. в схеме, которую Вы предлагаете реализовать, не укладываются входы DSP на последней стадии сложения если перенос с достаточной разрядностью для суммирования в DSP всего один, то такая схема не реализуемая
  7. по условию задачи и по обсуждению выше, в цифровой обработке допускается использовать БИХ фильтры второго порядка (только если нет высоких требований к точности коэффициентов и к точности вычислений) я отфильтровал сигналы из hilda_50kHz_rpm_high и hilda_50kHz_rpm_low в ВЧ фильтре с коэффициентами b = [1,-1]; a = [1,-0.9375]; с такой АЧХ: и затем отфильтровал в группе полосовых фильтров на разных частотах с такими коэффициентами: АЧХ фильтров: от результата фильтрации в группе полосовых фильтров взял по модулю и накопил процесс такого накопления в hilda_50kHz_rpm_high: результат такого накопления: результат накопления для hilda_50kHz_rpm_low: в сигнале hilda_50kHz_rpm_high больше всего накопил фильтр на частоте 3900. в сигнале hilda_50kHz_rpm_low больше всего накопил фильтр на частоте 400. если при других скоростях будут максимумы в фильтрах на соответствующих частотах, то таким образом можно будет приближенно вычислять скорость по частоте. если потребуется больше точность вычисления и будет позволять скорость вычислений, можно увеличить количество фильтров. больше фильтров - больше точность вычисления. если указанные коэффициенты не подходят по точности вычислений, то можно попробовать спроектировать фильтры попроще, но качество обработки может при этом снизиться.
  8. Здравствуйте! если рандомные "прострелы" длятся один такт и полоса полезного много меньше частоты дискретизации, то для очищения от таких рандомных "прострелов" может помочь медианный фильтр. я скачал файл hilda_50kHz_rpm_high.txt и попробывал применить медианный фильтр. оригинальный сигнал: результат обработки медианным фильтром третьего порядка: результат обработки медианным фильтром пятого порядка: по результатам обработки видно, что тут можно рассмотреть применение медианного фильтра до фильтрации ФНЧ и ФВЧ. к сожалению, в файле частота дискретизации уже пониженная до 50кГц. возможно, что применение медианного фильтра до понижения частоты лучше уберет рандомные "прострелы" и при этом меньше исказит полезный сигнал.
  9. если я правильно понимаю, по описанию работы шины генератор импульсов, счетчик и дешифратор (левый верхний в схеме) не могут переключаться во время передачи, так как иначе следующее в цепочке переключений устройство будет считать, что оно может передавать данные по общей шине. Как подобные коллизии решаются, если обратной связи устройств с генератором импульсов по схеме нет? что помешает счетчику переключиться во время передачи данных?
  10. представим, что у нас есть детализированная спектрограмма на длинном интервале времени, например вот такая 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); разглядеть детально небольшую область в таком масштабе действительно проблематично, поэтому выбираем Zoom In, Vertical Zoom любой участок спектрограммы приближается в пару кликов, все изменения графика вдоль оси времени - отлично отслеживаются для плавного скольжения по спектрограмме вдоль оси времени выбираем Pan, Vertical Pan довольно удобный способ исследования спектрограмм
  11. если уже есть вычисление данных и построение с помощью surf, например такое t = 0:1e-2:1; z = sin(2*pi*t).*cos(pi*t.'); surf(z); то получить двухмерный вариант можно например командой view(0,90), вот так t = 0:1e-2:1; z = sin(2*pi*t).*cos(pi*t.'); surf(z); view(0,90); таким образом будет выведен двухмерный вариант графика интенсивности, но так же сохранится возможность посмотреть и трехмерный график, с помощью кнопки Rotate 3D
  12. для проверки преобразования частоты с коэффициентом 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)); результаты моделирования такого фильтра получились похожими в самом Filter Designer были выбраны такие настройки
  14. пользуюсь сгенеренными в симулинке HDL файлами, впечатления хорошие очень похоже, что для Вашего случая подойдет блок FIR Rate Conversion HDL Optimized для проверки блока с дискретизацией 7/5 я использовал такие параметры результаты моделирования с входной дискретизацией как в Вашем случае (входной отсчет раз в два такта)
  15. функция importhdl обещает Import Verilog code and generate Simulink model Introduced in R2018b
  16. один из возможных способов расчета подобных выражений рассмотрим на примере суммы в числителе выражения посмотрим, как быстро сходится сумма для одного значения 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));
  17. эти вычисления - подстановка в передаточную функцию 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);
  18. в 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';
  19. ИХ двух фильтров каскадом - свертка ИХ отдельных фильтров. по рассчитанной ИХ можно построить АЧХ и ФЧХ для КИХ фильтров в матлабе можно рассчитать например так % 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); ИХ суммы/разности сигналов с двух фильтров - это сумма/разность ИХ отдельных фильтров. По рассчитанной ИХ можно построить АЧХ и ФЧХ
  20. вот пример, аналогичный тому что я пиводил ранее, но с учетом подсказки выше в обсуждении, о том что выгодно считать симметричную характеристику. % 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 - ИХ задержки на половину порядка рассчитанного фильтра. половина коэффициентов действительно получается нулевыми я не нашел, как в 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);
  21. для расчета активной и реактивной мощностей на известных гармониках, можно пойти и другим путем построим модель сигнала. Сначала сгенерим на нужной тактовой частоте нужную основную гармонику. Смоделируем вторую с заданной разностью фаз. (здесь я беру медленно линейно возрастающую до 60°). эти два сигнала будут имитировать принятые напряжение и ток добавим дополнительных кратных гармоник, добавим шум. в обработке сигнала перенесем с заданной известной частоты в ноль, при этом получится комплексный сигнал результат продецимируем, например в CIC фильтре по отфильтрованному продецимированному сигналу, найдем произведение напряжения и комплексно сопряженного с током накопим значения для очищения от шума (можно с помощью еще одно CIC фильтра) эта комплексная величина как-раз пропорциональна мощности гармоники: действительная часть пропорциональна активной составляющей, мнимая часть пропорциональна реактивной составляющей. теперь можем рассчитать значение коэффициента, на который нужно домножить активную мощность, чтобы получить реактивную делим чтобы убедиться, что алгоритм расчета устойчив к изменению входной частоты, изменим входную частоту на 58 Гц заметно изменился только сигнал после первого дециматора, частота ушла с центра на величину внесенной отстройки. в дальнейшем в алгоритме учтено что частота может быть отстроена, а начальная фаза неизвестна, на решение влияет именно разность фаз, поэтому дальше расчеты не изменились. такой алгоритм вычисляет отношение мощностей реактивной к активной только для одной гармоники, и поэтому дает хорошую оценку, только если большая часть мощности находится на основной гармонике. Если известно, что значимая для расчетов часть мощности находится на определенных кратных гармониках, то можно повторить этот алгоритм для этих определенных картых гармоник отдельными каналами, изменив только частоту переноса в начале. в алгоритме использованы простые по реализации фильтры и вычисления, для облегчения необходимой производительности. часть вычислений лежит за децимацией, что снижает требования к производительности в этой части вычислений.
  22. фильтр гильберта в Матлабе подходит для указанной задачи, но с некоторыми ограничениями. в задаче указано, что максимум порядок фильтра - несколько десятков, поэтому расчитаем фильтр гильберта для предельного случая, когда порядок фильтра равен 90, больше мы точно брать не можем. частота дискретизации 8кГц. Полоса фильтра, с учетом порядка фильтра, получилась от 50 Гц до 2кГц. для проверки расчитаных коэффициентов сравним передаточную функцию фильтра с передаточной функцией просто задержанного сигнала посмотрим ИХ построим АЧХ hfvt = fvtool(b0,1,b,1); hfvt.Fs = 8e3; hfvt.Analysis = 'magnitude'; построим ФЧХ hfvt = fvtool(b0,1,b,1); hfvt.Fs = 8e3; hfvt.Analysis = 'phase'; hfvt.PhaseUnits = 'degrees'; по результатам видно, что разность в 90° держится в необходимой полосе частот порядок фильтра равный 90 был взят как предельный возможный, при уменьшении порядка характеристики фильтра будут ухудшаться: сдвиг фазы все еще будет 90°, но полоса будет уменьшаться. для изучения влияния порядка фильтра на характеристики, сделаем аналогичные расчеты для более реального порядка фильтра, равного 30
  23. для вычисленных Вами квадратур формула atan2(Q,I) применяется и дает возможность увидеть исходную фазовую модуляцию я применял формулу к такому участку вычисленных Вами квадратур plot(t,I_signal,t,Q_signal); закон изменения фазы принятого сигнала plot(atan2(I_signal,Q_signal)); фаза вычисляется в пределах от минус пи до пи, для того, чтобы увидеть закон изменения фазы без скачков фазы при пересечении границ, можно воспользоваться функцией unwrap plot(unwrap(atan2(I_signal,Q_signal))); закон изменения фазы - порабола, как и положено для ЛЧМ сигнала. Также убедимся, что закон изменения частоты линейный plot(diff(unwrap(atan2(I_signal,Q_signal)))); по результатам расчетов видно, как именно использованный ранее в расчетах НЧ фильтр искажает закон изменения частоты сигнала при таких вычислениях
  24. надо проверять, и в данном случае можно разобраться на примере вашего фильтра тестбенч похож на предыдущий, но упрощен входной сигнал сначала посмотрим вариант с такой конструкцией always @(posedge clk) begin fifo_sum <= fifo_sum + datain - FIFO[k]; end always @(negedge clk) begin FIFO[k] <= datain; k <= k + 1; end фильтр работает не так как задумывалось накопление в фильтре идет, но размер окна равен 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 все работает так как предполагалось изначально счетчик (k) расчитывается в @(posedge clk) и используется в @(posedge clk), кольцевой буфер дает задержку ровно на 16 тактов, как и положено, размер окна равен 16