如何在 Arduino 中进行快速互相关（实时即时信号识别器的一部分）
How to do Fast Cross-Correlation in Arduino (Part of Real Time Instant Signal Recognizer)
我很难编写代码来对信号（在引脚 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;
}
感谢帮助的人
我会尽量做到通用 - 解决方案仍然适用于您的问题。
注意：我希望优化方程式，而不是代码。我不明白代码与等式的关系。
假设您只需要找到信号的最后 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 不会更改，请将其完全从公式中删除。提取零只是浪费时间，没有任何好处。
另一个技巧，如果应用程序允许，将降低采样率。这样，样本（以及因此中断）变得更慢，处理器有更多时间来完成这项工作。
我很难编写代码来对信号（在引脚 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;
}
感谢帮助的人
我会尽量做到通用 - 解决方案仍然适用于您的问题。
注意：我希望优化方程式，而不是代码。我不明白代码与等式的关系。
假设您只需要找到信号的最后 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 不会更改，请将其完全从公式中删除。提取零只是浪费时间，没有任何好处。
另一个技巧，如果应用程序允许，将降低采样率。这样，样本（以及因此中断）变得更慢，处理器有更多时间来完成这项工作。