从手机的加速度计检测周期性数据

Question

从手机的加速度计检测周期性数据

androidalgorithmsignal-processingaccelerometer

4

我正在开发一款Android应用程序，需要检测用户的上下文（是否在最小行驶或步行）。

我正在使用加速度计和所有轴的总和来检测加速度矢量。在步行时，它的工作效果相当不错，我可以看到一些周期性的值。但是我需要以编程方式检测这些周期。

请问是否有任何数学函数可以检测一组值中的周期？我听说傅里叶变换可用于此，但我真的不知道如何实现它。它看起来非常复杂：）

请帮帮我

- simekadam

3个回答

7

还有一种使用FFT计算您的数据自相关性的方法，将复杂度从O(n^2)降低到O(n log n)。基本思想是将周期样本数据转换为FFT，然后通过将每个FFT系数乘以其复共轭来计算功率谱，然后对功率谱进行反FFT。你可以很容易地找到现成的代码来计算功率谱。例如，看看Moonblink android library。该库包含FFTPACK的JAVA翻译（一个好的FFT库），并且还有一些用于计算功率谱的DSP类。我成功使用的自相关方法是McLeod Pitch Method（MPM），其java源代码可在here中找到。我编辑了类McLeodPitchMethod中的一个方法，使其能够使用FFT优化的自相关算法计算音高：

private void normalizedSquareDifference(final double[] data) {
    int n = data.length;
    // zero-pad the data so we get a number of autocorrelation function (acf)
    // coefficients equal to the window size
    double[] fft = new double[2*n];
        for(int k=0; k < n; k++){
        fft[k] = data[k];
    }
    transformer.ft(fft);
    // the output of fft is 2n, symmetric complex
    // multiply first n outputs by their complex conjugates
    // to compute the power spectrum
    double[] acf = new double[n];
    acf[0] = fft[0]*fft[0]/(2*n);
    for(int k=1; k <= n-1; k++){
        acf[k] = (fft[2*k-1]*fft[2*k-1] + fft[2*k]*fft[2*k])/(2*n);
    }
    // inverse transform
    transformerEven.bt(acf);
    // the output of the ifft is symmetric real
    // first n coefficients are positive lag acf coefficients
    // now acf contains acf coefficients
    double[] divisorM = new double[n];
    for (int tau = 0; tau < n; tau++) {
        // subtract the first and last squared values from the previous divisor to get the new one;
        double m = tau == 0 ? 2*acf[0] : divisorM[tau-1] - data[n-tau]*data[n-tau] - data[tau-1]*data[tau-1];
        divisorM[tau] = m;
        nsdf[tau] = 2*acf[tau]/m;
    }
}

其中transformer是来自Java FFTPACK翻译的FFTTransformer类的私有实例，而transformerEven是FFTTransformer_Even类的私有实例。使用您的数据调用McLeodPitchMethod.getPitch()将给出非常有效的频率估计。

- willpett

我的一个问题是：如果您有一个简单的FFT可用，那么在检测最强周期分量时使用它来计算自相关的目的是什么？直接从FFT中获取最强的成分并跳过其余部分不是可以吗？就我所理解的，自相关对于这个任务只是一种不需要理解和实现FFT的方法，因此如果已经有FFT可用，我不再看到它的意义。 - LiKao

1

这是一种策略，确实如此。然而，在FFT频谱中识别“最强组件”并不总是直截了当的。具有最大幅度的分量可能不代表数据的基本频率，或者更糟糕的是，基本频率可能根本不存在于频谱中。虽然自相关也不能直接测量基本频率，但在面对这些问题时，它是最好的估计方法。如果您怀疑这些问题不适用于您的数据，则可能认为自相关是不必要的。 - willpett

0

这是一个使用libgdx的FFT类来计算安卓傅里叶变换的示例：

    package com.spec.example;
import android.app.Activity;
import android.os.Bundle;
import com.badlogic.gdx.audio.analysis.FFT;
import java.lang.String;
import android.util.FloatMath;
import android.widget.TextView;

public class spectrogram extends Activity {
    /** Called when the activity is first created. */
    float[] array = {1, 6, 1, 4, 5, 0, 8, 7, 8, 6, 1,0, 5 ,6, 1,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    float[] array_hat,res=new float[array.length/2];
    float[] fft_cpx,tmpr,tmpi;
    float[] mod_spec =new float[array.length/2];
    float[] real_mod = new float[array.length];
    float[] imag_mod = new float[array.length];
    double[] real = new double[array.length];
    double[] imag= new double[array.length];
    double[] mag = new double[array.length];
    double[] phase = new double[array.length];
    int n;
    float tmp_val;
    String strings;
    FFT fft = new FFT(32, 8000);
    @Override
    public void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        TextView tv = new TextView(this);

       fft.forward(array);
       fft_cpx=fft.getSpectrum();
       tmpi = fft.getImaginaryPart();
       tmpr = fft.getRealPart();
      for(int i=0;i<array.length;i++)
       {
           real[i] = (double) tmpr[i];
           imag[i] = (double) tmpi[i];
           mag[i] = Math.sqrt((real[i]*real[i]) + (imag[i]*imag[i]));
           phase[i]=Math.atan2(imag[i],real[i]);

           /****Reconstruction****/        
           real_mod[i] = (float) (mag[i] * Math.cos(phase[i]));
           imag_mod[i] = (float) (mag[i] * Math.sin(phase[i]));

       }

       fft.inverse(real_mod,imag_mod,res);
   }
}

更多信息请参见：http://www.digiphd.com/android-java-reconstruction-fast-fourier-transform-real-signal-libgdx-fft/

- Manu

1

请注意，这只是使用FFT的简单示例，无法描述FFT本身的实际实现/算法。这更像是libgdx的单元测试，如果没有深入了解傅里叶变换，它就毫无价值。正如页面上所说：“请随意使用此代码，但在使用之前，请至少学习一下它的工作原理。” - LiKao

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- LiKao · Accepted Answer

检测数据周期的最简单方法是自相关。这也相当容易实现。要获取在i处的自相关性，只需将数据中的每个数据点与移位i的每个数据点相乘即可。以下是一些伪代码：

for i = 0 to length( data ) do
  autocorrel[ i ] = 0;
  for j = 0 to length( data ) do
     autocorrel[ i ] += data( j ) * data( ( j + i ) mod length( data ) )
  done
done

这将给您一个值数组。最高的“周期性”在具有最高值的索引处。这样，您可以提取任何周期部分（通常不止一个）。

此外，我建议您不要尝试在应用程序中实现自己的FFT。虽然这个算法非常适合学习，但是有很多错误可能会出现，并且很难进行测试，而且您的实现也很可能比已经可用的实现慢得多。如果您的系统允许，我建议您使用FFTW，在FFT实现方面无论如何都无法超越。

编辑：解释为什么这即使在没有完全重复的值上也有效。

通常而完全正确的计算自相关的方法是从数据中减去平均值。比如说你有[1, 2, 1.2, 1.8]，然后你可以从每个样本中提取1.5，留下[-.5, .5, -.3, .3]。现在，如果你将其在偏移量为零的位置相乘，负数将与负数相乘，正数将与正数相乘，得到(-.5)^2 + (.5)^2 + (-.3)^2 + (.3)^2=.68。在偏移量为1的位置，负数将与正数相乘，得到(-.5)*(.5) + (.5)*(-.3) + (-.3)*(.3) + (.3)*(-.5)=-.64。在偏移量为2的位置，负数再次与负数相乘，正数与正数相乘。在偏移量为3的位置，类似于偏移量为1的情况。你会发现，在偏移量为0和2（周期）时，得到的是正值，在1和4处得到的是负值。

现在仅需检测周期，无需减去平均值。如果您只是将样本保留原样，则每次添加时都会添加平方平均值。由于对于计算的每个系数都会添加相同的值，因此与首先减去平均值相比，比较将产生相同的结果。最坏的情况是数据类型可能会溢出（如果使用某种整数类型），或者当值开始变得很大时可能会出现舍入误差（如果使用float，通常不是问题）。如果发生这种情况，请先减去平均值，然后尝试是否可以获得更好的结果。

使用自相关与某种快速傅里叶变换相比，最大的缺点是速度。自相关需要O(n^2)，而FFT仅需要O(n log(n))。如果您需要经常计算非常长序列的周期，则自相关可能不适用于您的情况。

如果您想了解傅里叶变换的工作原理，以及关于实部、虚部、幅度和相位的所有内容（例如Manu发布的代码），建议您查看这本书。

编辑2：在大多数情况下，数据既不完全周期性，也不完全混沌和非周期性。通常，您的数据将由多个周期组成，强度各异。周期是一种时间差，您可以通过它来移动数据，使其类似于自身。自相关计算数据的相似度，如果您将其移位某个量。因此，它为您提供了所有可能周期的强度。这意味着，没有“重复值索引”，因为当数据完全周期时，所有索引都将重复。具有最强值的索引给出了时间偏移量，而不是数据中的索引。为了理解这一点，重要的是要理解时间序列如何被认为是由完全周期函数（正弦基函数）的总和组成的。

如果您需要检测长时间序列，则通常最好在数据上滑动一个窗口，仅检查较小数据框的周期。但是，您必须意识到，您的窗口将向您的数据添加额外的周期，您必须注意这一点。

更多信息请参见我在上次编辑中发布的链接。