Получение «пиков» из массива значений — измерение оборотов с помощью акселерометра

Question

Получение «пиков» из массива значений — измерение оборотов с помощью акселерометра

arduino-nano

Для проекта, в котором я собираюсь использовать Arduino и акселерометр для определения оборотов дизельного/бензинового двигателя, я намерен использовать следующее решение. Я приклеил акселерометр mma8453 к неодовому магниту, который будет размещен на металлической части двигателя. Используя Arduino, я читаю значения с акселерометра и вычисляю одно конечное значение, которое равно sqrt (acc.x^2 + acc.y^2 + acc.z^2), а также получаю миллисекунды для каждого значения. Значения выглядят следующим образом:

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

Моя проблема в том, что я понятия не имею, как написать функцию для этого.

Любая помощь действительно приветствуется

Спасибо

269503  310     
269504  303     
269506  334     
269507  326     
269509  331     
269510  330     
269512  380     
269513  377     
269514  392     
269516  400     
269517  339     
269519  393     
269520  368     
269521  353     
269523  346     
269524  316     
269527  330     
269528  281     
269530  196     
269531  195     
269532  162     
269534  119     
269535  140     
269537  132     
269538  90      
269539  99      
269541  122     
269542  168     
269544  224     
269545  219     
269546  211     
269548  268     
269549  301     
269551  334     
269552  332     
269554  316     
269555  329     
269556  379     
269558  363     
269559  381     
269561  388     
269562  403     
269564  400     
269565  420 49  
269566  388     
269569  360     
269570  346     
269572  259     
269573  276     
269574  284     
269576  242     
269577  207     
269579  192     
269580  145     
269582  119     
269583  113     
269584  45      
269586  125     
269587  142     
269589  123     
269590  205     
269591  243     
269593  304     
269594  308     
269596  305     
269597  298     
269599  321     
269600  340     
269601  360     
269603  380     
269604  384     
269606  370     
269607  372     
269608  413     
269611  441 46  
269612  395     
269614  355     
269615  332     
269617  290     
269618  317     
269619  286     
269621  227     
269622  190     
269624  172     
269625  199     
269627  146     
269628  154     
269629  177     
269631  176     
269632  156     
269634  136     
269635  174     
269636  212     
269638  260     
269639  318     
269641  291     
269642  321     
269644  304     
269645  315     
269646  343     
269648  339     
269649  360     
269651  393     
269652  390     
269654  305     
269656  405     
269657  458 46  
269659  358     
269660  372     
269662  339     
269663  277     
269664  298     
269666  241     
269667  203     
269669  208     
269670  185     
269672  130     
269673  165     
269674  105     
269676  87      
269677  109     
269679  128     
269680  176     
269681  229     
269683  226     
269684  271     
269686  343     
269687  307     
269688  320     
269690  318     
269691  322     
269693  322     
269694  369     
269697  363     
269698  413     
269699  399     
269701  388     
269702  426 45  
269704  395     
269705  380     
269707  342     
269708  311     
269709  309     
269711  315     
269712  221     
269714  227     
269715  198     
269716  113     
269718  111     
269719  88      
269721  44      
269722  69      
269723  107     
269725  158     
269726  200     
269728  232     
269729  247     
269731  313     
269732  338     
269733  348     
269735  366     
269736  339     
269739  348     
269740  377     
269741  378     
269743  385     
269744  394     
269746  403     
269747  447     
269749  454 47  
269750  398     
269751  380     
269753  302     
269754  245     
269756  299     
269757  310     
269758  238     
269760  206     
269761  180     
269763  108     
269764  81      
269766  83      
269767  11      
269768  119     
269770  106     
269771  169     
269773  245     
269774  270     
269775  311     
269777  338     
269778  341     
269780  332     
269782  337     
269783  337     
269785  380     
269786  387     
269788  382     
269789  373     
269791  398     
269792  414 43  
269793  408     
269795  387     
269796  374     
269798  285     
269799  290     
269801  326     
269802  266     
269803  211     
269805  221     
269806  191     
269808  120     
269809  175     
269810  63      
269812  147     
269813  59      
269815  120     
269816  183     
269817  196     
269819  252     
269820  301     
269822  334     
269824  325     
269826  334     
269827  310     
269828  351     
269830  354     
269831  358     
269833  373     
269834  379     
269835  332     
269837  386     
269838  404 46  
269840  383     
269841  385     
269843  332     
269844  266     
269845  301     
269847  282     
269848  188     
269850  199     
269851  182     
269852  128     
269854  163     
269855  119     
269857  120     
269858  140     
269860  120     
269861  129     
269862  212     
269864  230     
269865  241     
269868  290     
269869  304     
269871  340     
269872  350     
269873  327     
269875  323     
269876  351     
269878  369     
269879  387     
269880  379     
269882  380     
269883  421 45  
269885  373     
269886  364     
269888  360     
269889  305     
269890  302     
269892  266     
269893  249     
269895  202     
269896  229     
269897  136     
269899  122     
269900  128     
269902  49      
269903  98      
269904  107     
269906  158     
269907  211     
269910  235     
269911  225     
269913  269     
269914  314     
269915  341     
269917  335     
269918  325     
269920  348     
269921  410     
269923  375     
269924  404     
269925  419     
269927  425     
269928  442 45  
269930  436     
269931  400     
269932  393     
269934  326     
269935  259     
269937  307     
269938  265     
269940  220     
269941  199     
269942  183     
269944  86      
269945  83      
269947  98      
269948  23      
269949  65      
269952  124     
269953  152     
269955  221     
269956  249     
269957  287     
269959  350     
269960  331     
269962  323     
269963  341     
269964  338     
269966  389     
269967  388     
269969  380     
269970  370     
269972  358     
269973  415 45  
269974  403     
269976  377     
269977  381     
269979  338     
269980  328     
269982  284     
269983  246     
269984  197     
269986  178     
269987  230     
269989  162     
269990  171     
269991  141     
269993  258     
269995  138     
269997  111     
269998  163     
270000  212     
270001  192     
270002  304     
270004  287     
270005  294     
270007  349     
270008  309     
270009  334     
270011  332     
270012  364     
270014  361     
270015  408     
270017  336     
270018  421     
270019  432 46  
270021  361     
270022  361     
270024  342     
270025  271     
270027  288     
270028  310     
270029  235     
270031  218     
270032  166     
270034  136     
270035  146     
270036  102     
270039  78      
270040  99      
270042  122     
270043  105     
270044  181     
270046  222     
270047  293     
270049  299     
270050  299     
270052  318     
270053  349     
270054  344     
270056  349     
270057  382     
270059  381     
270060  412     
270062  380     
270063  397     
270064  425 45  
270066  381     
270067  367     
270069  332     
270070  272     
270071  291     
270073  296     
270074  178     
270076  177     
270077  190     
270080  123     
270081  146     
270082  107     
270084  85      
270085  108     
270087  117     
270088  148     
270089  234     
270091  223     
270092  238     
270094  283     
270095  301     
270097  332     
270098  331

@lucian_v, 👍-1

Обсуждение

(Это домашнее задание или вы действительно что-то проектируете?) Вы упускаете такие вещи, как почему вы используете магнит? Где вы размещаете часть NXP? Это оставляет нам возможность догадываться о том, что происходит. Рассчитав на основе предположений, я бы сказал, что ваш двигатель работает со скоростью около 1304,3 об/мин., @st2000

На самом деле это смесь домашней работы и дизайна. Акселерометр крепится на магните и вместе с ним крепится к двигателю. Собственно измерьте вибрации двигателя. Когда есть детонация, то у меня высокое значение. Я не знаю, что вы имеете в виду под NXP., @lucian_v

NXP (также известная как Motorola, Freescale и, возможно, Qualcomm) является OEM-производителем указанного вами чипа., @st2000

3 ответа

@*tttapa* · Answer 1 · 2019-06-22 17:10-00:00

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

Вот пример, который работает с вашими данными. Он написан на Python для удобства тестирования и построения графиков, но транскодирование его на C++ должно быть тривиальным.

График исходных данных, отфильтрованных данных с использованием EMA и рассчитанного RPM

Я использовал Фильтр низких частот Exponential Moving Average, чтобы избавиться от шума. В зависимости от задействованных частот вы можете изменить коэффициент сдвига. Чем выше этот коэффициент, тем ниже частота среза и тем больше фазовый сдвиг.

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

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

Вы можете добавить еще одно EMA, чтобы сгладить измерения RPM.

class ZeroCrossRPMCounter:
    def __init__(self,
                 max_rpm: float):
        self.min_distance = int(round(60_000 / max_rpm))

        self.first_crossing = True
        self.time_of_last_crossing = 0
        self.rpm = float('nan')
        self.previous_value = 0

    def __call__(self, time_ms: int, value: int):
        new_crossing = self.previous_value >= 0 and value < 0
        if new_crossing:
            if self.first_crossing:
                self.time_of_last_crossing = time_ms
                self.first_crossing = False
            elif (time_ms - self.time_of_last_crossing) >= self.min_distance:
                self.rpm = 60_000 / (time_ms - self.time_of_last_crossing)
                self.time_of_last_crossing = time_ms
            else:
                new_crossing = False
        self.previous_value = value
        return new_crossing


class EMA_shift_round:
    def __init__(self, shiftfac: int, initial: int = 0):
        self.shiftfac = shiftfac
        self.fixedPointOneHalf = 1 << (shiftfac - 1)
        self.z = initial << shiftfac

    def __call__(self, x: int):
        self.z += x
        shifted = (self.z + self.fixedPointOneHalf) >> self.shiftfac
        self.z -= shifted
        return shifted


if __name__ == '__main__':
    import numpy as np
    from matplotlib import pyplot as plt

    # Load the data
    data = np.genfromtxt('rpm-ctr-data.csv', delimiter='\t', dtype=np.int)
    time, signal = data[:, 0], data[:, 1]
    t_start, t_end = time[0] * 1e-3, time[-1] * 1e-3

    # Find the DC offset
    offset = int(round(np.mean(signal)))
    print(offset)

    # Initialize the filter and RPM counter
    lp_filter = EMA_shift_round(3)
    counter = ZeroCrossRPMCounter(2000)

    # Iterate over the data, filter it, and count the RPM
    peaks = []
    rpms = []
    filtered = np.zeros(len(time))
    for n in range(len(time)):
        filtered[n] = lp_filter(signal[n] - offset)
        new_crossing = counter(time[n], filtered[n])
        if new_crossing:
            peaks.append(time[n] * 1e-3)
            rpms.append(counter.rpm)

    # Plot the results
    plt.figure(figsize=(10, 6))
    plt.subplot(3, 1, 1)
    plt.title('Original signal')
    plt.plot(time * 1e-3, signal)
    for peak in peaks:
        plt.axvline(peak, color='red')
    plt.xlim((t_start, t_end))
    plt.subplot(3, 1, 2)
    plt.title('Filtered signal')
    plt.plot(time * 1e-3, filtered)
    for peak in peaks:
        plt.axvline(peak, color='red')
    plt.axhline(0, color='black', linestyle=':')
    plt.xlim((t_start, t_end))
    plt.subplot(3, 1, 3)
    plt.title('RPM measurements')
    plt.plot(peaks, rpms, '-o')
    plt.axhline(np.mean(rpms[1:]), color='red', linestyle=':')
    plt.xlim((t_start, t_end))
    plt.xlabel("Time [s]")
    plt.tight_layout()
    plt.show()

@*st2000* · Answer 2 · 2019-06-22 20:11-00:00

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

Тем не менее, использование БПФ (быстрых преобразований Фурье) — это альтернативный, но принципиально похожий подход к обнаружению пересечений нуля для нахождения основные частоты и, следовательно, обороты двигателя.

Это " БПФ на Arduino" описывает, как реализовать быстрое преобразование Фурье на Arduino. Он использует библиотеки БПФ, размещенные здесь, на github.

@*Edgar Bonet* · Answer 3 · 2019-06-22 20:22-00:00

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

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

const int LOW_THRESHOLD = 220;
const int HIGH_THRESHOLD = 300;

/* Возвращает true, если обнаружено превышение высокого порога. */
void crossing_detector(int x)
{
    static enum { BELOW, ABOVE } state;
    if (state == BELOW && x >= HIGH_THRESHOLD) {
        state = ABOVE;
        return true;
    }
    if (state == ABOVE && x < LOW_THRESHOLD) {
        state = BELOW;
    }
    return false;
}

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