Как сделать быструю взаимную корреляцию в Arduino (часть распознавателя мгновенных сигналов в реальном времени)

Мне трудно написать код, который отбирает сигнал (на выводе A0) и выполняет мгновенную перекрестную корреляцию в реальном времени с другим известным сигналом (который сохраняется во флэш-памяти Arduino Uno).

Моя проблема в том, что мой код (мое уравнение) включает в себя множество циклов, из-за чего Arduino сложно делать это в реальном времени и сразу после выборки сигнала (используя ISR (TIMER1_COMPA_vect)).

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

Уравнение взаимной корреляции: введите здесь описание изображения

для протокола, вот что я уже сделал:

#define NUM_OF_SAMPLES 64
#define NUM_OF_ETALONS 6

const double* const Etalon_Signals[NUM_OF_ETALONS] PROGMEM = {sinWave1, sinWave2, sinWave3, rectWave4, rectWave5, triWave6};
const double Etalon_Signals_Mean[NUM_OF_ETALONS] PROGMEM = {sinWaveMean1, sinWaveMean2, sinWaveMean3, rectWaveMean4, rectWaveMean5, triWaveMean6};
const double Etalon_Signals_Variance[NUM_OF_ETALONS] PROGMEM = {sinWaveVariance1, sinWaveVariance2, sinWaveVariance3, rectWaveVariance4, rectWaveVariance5, triWaveVariance6};
//-------------- SAMPLED SIGNAL ---------------------
volatile double SignalSamples[NUM_OF_SAMPLES];
volatile int idx = 0;
//-------------- CORRELATION VARS -------------------
// Correlations[NUM_OF_ETALONS] -         array of 6 values which calculate the correlation for particular 'd' (delay) value.
// MaxCorrelationValues[NUM_OF_ETALONS] - array which holds the maximum value of a praticular Etalon signal,
//                                        which I use later to decide which wave is most simmilar to the sampled wave.
//                                        MaxCorrelationValues[k] = the maximum correlation value out of all the 'd' iterations, 
//                                        of Etalon wave number k (look at the formula for better understanding).
// MostSimillarCorrVal -                  Holds the Maximum value of correlation out of 6 Etalons (every Etalon has a maximum correlation point out of d correlations)
//                                        MostSimillarCorrVal will hold the Absulote Maximum out of the 6 Etalons.
// SignalSamplesVar -                     Variance of the Sampled Signal -> Eq: (sum(y[i] - my))
// SignalSamplesMean -                    Mean of the Sampled Signal -> Eq: (my)
// counter -                              used for smart Circular scanning over the ring buffer "SignalSamples".
// maxCorrWave -                          Holds the 2^index of the most simmilar wave out of the etalons.
// SignalSamplesVarTemp -                 a particular substraction of the signal with its mean value -> Eq: (numerator) -> (y(i-d)-my).
volatile double Correlations[NUM_OF_ETALONS] = {0}, MaxCorrelationValues[NUM_OF_ETALONS] = {0}, MostSimillarCorrVal = 0;
volatile double SignalSamplesVar = 0, SignalSamplesMean = 0;
volatile int counter = 0, maxCorrWave = 0;
double EtalonVarianceVal[NUM_OF_ETALONS], SignalSamplesVarTemp;


void setup()
{
  SetTimer1CTCforFreq(SAMPLING_FREQUENCY);
  DDRB |= WAVE_1 | WAVE_2 | WAVE_3 | WAVE_4 | WAVE_5 | WAVE_6;
  pinMode(A0, INPUT);
  for(int i = 0; i<NUM_OF_ETALONS;i++)
    EtalonVarianceVal[i] = pgm_read_float(&Etalon_Signals_Variance[i]);
}

ISR(TIMER1_COMPA_vect)
{
  SignalSamplesVar -= (SignalSamples[idx] - SignalSamplesMean)*(SignalSamples[idx] - SignalSamplesMean);
  SignalSamplesMean = SignalSamplesMean * NUM_OF_SAMPLES - SignalSamples[idx];
  SignalSamples[idx] = analogRead(ANALOG_INPUT_PIN);                        // sampling
  SignalSamplesMean = (SignalSamplesMean + SignalSamples[idx])/NUM_OF_SAMPLES;
  SignalSamplesVar +=  (SignalSamples[idx] - SignalSamplesMean)*(SignalSamples[idx] - SignalSamplesMean);
  // SignalSamplesMean = Mean(SignalSamples, NUM_OF_SAMPLES);                  // Mean of current state calculations.
  // SignalSamplesVar = VarianceSquared(SignalSamples, NUM_OF_SAMPLES);        // Variance of current state calculations.
  for(int delay = -NUM_OF_SAMPLES; delay < NUM_OF_SAMPLES; delay++)         // "shifting" array loop.
  {
    counter = 0;                                                            // counter - an index which runs on SignalSamples array.
    while(counter < NUM_OF_SAMPLES)                                         // loop running on the SignalSamples in circle.
    {
      // SignalSamplesVarTemp = Signal's "Variance" semi-calculation.
      SignalSamplesVarTemp = SignalSamples[(idx + counter - delay)%NUM_OF_SAMPLES] - SignalSamplesMean;
      // Calculating all Cross Correlations with all the six Etalons.
      for(int k = 0; k < NUM_OF_ETALONS; k++)
        Correlations[k] += pgm_read_float(&Etalon_Signals[k][counter]) * SignalSamplesVarTemp; // eMean = 0 always.
      counter++;
    }
    for(int k = 0; k < NUM_OF_ETALONS; k++) // updating all the Etalons.
    {
      Correlations[k] *= Correlations[k];
      Correlations[k] /= SignalSamplesVar * pgm_read_float(&EtalonVarianceVal[k]);        // The Correlation for particular delay.
      // MaxCorrelationValues[k] - Max value of the correlation with Etalon[k]. 
      // At the end we compare which etalon got the biggest correlation value (which is positive).
      MaxCorrelationValues[k] = MaxCorrelationValues[k] < Correlations[k] ? Correlations[k] : MaxCorrelationValues[k]; 
    }
  }
  for(int k = 0; k < NUM_OF_ETALONS; k++)
  {
    if(MaxCorrelationValues[k] > MostSimillarCorrVal)
    {
      MostSimillarCorrVal = Correlations[k];
      maxCorrWave = k;
    }
  }
  PORTB = 1 << maxCorrWave;
  idx = (idx+1)%NUM_OF_SAMPLES;
}

Спасибо тем, кто помогает


person Jhon Margalit    schedule 20.07.2020    source источник
comment
Насколько я знаю, выполнять тяжелые вычисления внутри ISR всегда плохая идея. В ISR вы просто сохраняете важные данные, а затем обрабатываете эти данные вне ISR. Как? Это зависит от выбранной архитектуры/дизайна.   -  person virolino    schedule 20.07.2020
comment
Если числа достаточно малы и нет риска переполнения, вы можете извлечь квадратный корень из произведения, и это сэкономит вам много вычислительной мощности и времени.   -  person virolino    schedule 20.07.2020
comment
Наверняка именно поэтому в коде нет функции sqrt(). Это все еще плохо работает. Я не знаю, как минимизировать количество циклов.   -  person Jhon Margalit    schedule 20.07.2020
comment
Проблема выполнения тяжелых вычислений внутри ISR все еще остается. В любом случае, дело в том, что я не вижу, как предоставленный вами код реализует приведенную выше формулу - возможно, потому, что корреляции - это не то, что я использую каждый день (читай: совсем не за последние 20 с лишним лет). Если вы немного объясните, что вы сделали, я мог бы быть более полезным.   -  person virolino    schedule 20.07.2020


Ответы (1)


Я постараюсь быть универсальным - решение по-прежнему применимо к вашей проблеме.

ПРИМЕЧАНИЕ. Я стараюсь оптимизировать уравнение, а не код. Я не понимаю, как код относится к уравнению.

Предположим, вам нужно найти только сумму последних N отсчетов сигнала.

Пусть образцы будут

x0, x1, ..., x(N-1)

Пусть первая сумма

S0 = x0 + x1 + ... + x(N-1)

Затем прибывает новый образец:

xN

Сумма

S1 = x1 + x2 + ... + x(N-1) + xN

Но эту сумму можно переписать как:

S1 = S0 - x0 + xN

Таким образом, вам нужно запомнить всего 3 числа:

  • предыдущая/текущая сумма (S0);
  • самое старое значение выборки (x0);
  • новейшее значение выборки (xN).

Поэтому вместо вычисления (N-2) сложений вы вычисляете только 2 (фактически только сложение и вычитание).


Если вам нужно вычислить сумму

x1 - xm, x2 - xm ...

вы должны заметить, что при расчете S1 у вас будет

S1 = S0 - (x0-xm) + (xN-xm)

что переводится как

S1 = S0 - x0 + xm

что избавляет вас от добавления и вычитания xm без необходимости.


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

Кроме sqrt() - я понятия не имею, как тут помочь.


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

На самом деле, первый S0 будет

S0 = x0

Другой оптимизацией будет вычисление

y(i-d) - my

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


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

person virolino    schedule 20.07.2020
comment
Вы правы, но пока ничего нового. самые первые строки в ISR делают именно то, что вы только что сказали. sqrt() можно удалить, так как я возводил в квадрат все уравнение, поэтому в sqrt() нет необходимости. моя проблема в том, что я не знаю, смогу ли я удалить первый цикл, который проходит через переменную «задержка» (в уравнении это d). - person Jhon Margalit; 20.07.2020
comment
Видеть? Наконец-то ты начал делать правильные вещи. Пожалуйста, отредактируйте вопрос и добавьте туда эту информацию. Кроме того, внимательно опишите: уравнения, которые вы собираетесь использовать, псевдокод того, что вы хотите сделать (без оптимизации), < b>оптимизации, которые вы имеете в виду, и реальный код, включая оптимизации. После этого мы можем проанализировать, что вы еще можете сделать, если что-то еще можно сделать. Если они не очевидны, также опишите, как вышеперечисленные вещи соотносятся друг с другом. Например, какая переменная в коде относится к mx в исходном уравнении. - person virolino; 20.07.2020
comment
Разве значения mx и my не меняются с каждым новым образцом? - person the busybee; 20.07.2020
comment
Может быть, они меняются, но я не знаю, как, может быть, ОП может помочь... - person virolino; 20.07.2020
comment
'''my''' меняется, '''mx''' известно (это среднее значение сигнала эталона, которое сохраняется во флэш-памяти (и в моем случае это тоже 0). Поскольку X[] является эталоном, он постоянен, поэтому mx не меняется. - person Jhon Margalit; 21.07.2020
comment
@JhonMargalit: Кстати, вы уверены, что у вашего процессора достаточно мощности, чтобы делать то, что вы хотите? - person virolino; 22.07.2020
comment
@JhonMargalit: как насчет my? какова его формула (если он непостоянен)? - person virolino; 22.07.2020
comment
@virolino my - это среднее значение входного сигнала. Формула такова: SUM(y)/LENGTH_OF_Y, просто среднее. но я не буду делать цикл над этим, потому что это драгоценное время, поэтому я делаю каждую ISR (прерывание). Я вычитаю самое старое значение y из my, а затем, когда я делаю выборку на этой итерации, я добавляю выборку к среднему. это намного быстрее, чем выполнение цикла по выборкам y на каждой итерации ISR (то есть O(n)). А насчет мощности процессора да, уверен, что можно. - person Jhon Margalit; 22.07.2020