Matlab:高效生成听觉尖峰电位

Question

我正在尝试在MATLAB中从采样信号生成一系列“听觉尖峰”，到目前为止我实现的方法很慢。太慢了。

尖峰是由http://patrec.cs.tu-dortmund.de/pubs/papers/Plinge2010-RNF.pdf的第II，B部分产生的:检测过零点，并对每对信号之间的(平方根压缩)正部分求和。这给出了每个峰值的高度。它们的位置是通过识别每对零交叉之间具有最大值的样本来找到的。

我想过使用accumarray(...)为此，我产生了生成一个矩阵的想法，其中每一列代表一对零交叉(一个尖峰)，每一行代表一个样本。然后，用相应的一对零交叉之间的1填充每一列。

当前的实现从实际的数据向量中填充这些列，这样我们以后就不必使用accumarray。

当前实施：

function out = audspike(data)
    % Find the indices of the zero crossings. Two types of zero crossing:
    %   * Exact, samples where data == 0
    %   * Change, where data(i) .* data(i+1) < 0; that is, data changes sign
    %   between two samples. In this implementation i+1 is returned as the
    %   index of the zero crossing.
    zExact = (data == 0);
    zChange = logical([0; data(1:end-1) .* data(2:end) < 0]);
    zeroInds = find(zExact | zChange);

    % Vector of the difference between each zero crossing index
    z=[zeroInds(1)-1; diff(zeroInds)];

    % Find the "number of zeros" it takes to move from the first sample to the 
    % a given zero crossing
    nzeros=cumsum(z);

    % If the first sample is positive, we cannot generate a spike for the first
    % pair of zero crossings as this represents part of the signal that is
    % negative; therefore, skip the first zero crossing and begin pairing from
    % the second
    if data(1) > 0
        nzeros = nzeros(2:2:end);
        nones = z(3:2:end)+1;
    else
        nzeros = nzeros(1:2:end);
        nones = z(2:2:end)+1;
    end

    % Allocate sparse array for result
    G = spalloc(length(data), length(nzeros), sum(nones));

    % Loop through pairs of zero crossings. Each pair gets a column in the
    % resultant matrix. The number of rows of the matrix is the number of 
    % samples. G(n, ii) ~= 0 indicates that sample n belongs to pair ii
    for ii = 1:min(length(nzeros), length(nones))
        sampleInd = nzeros(ii)+1:nzeros(ii)+nones(ii)-1;
        G(sampleInd, ii) = data(sampleInd);
    end

    % Sum the square root-compressed positive parts of signal between each zero
    % crossing
    height = sum(sqrt(G), 2);

    % Find the peak over average position
    [~, pos] = max(G, [], 2);

    out = zeros(size(data));
    out(pos) = height;
end

正如我所说的，这很慢，而且它一次只适用于一个通道。最慢的部分(毫不奇怪)是循环。如果我将矩阵G的分配更改为标准零(...)然后，慢的部分变成sum(...)，而不是稀疏数组和max(...)计算，原因很明显。

我怎样才能更有效率地做到这一点？我并不反对写一个MEX函数，如果需要的话。

Answer 1

我现在已经创建了一个MEX函数，它执行与上面的代码相同的任务，除了最后三行之外，这三行代码只需要很少的修改就可以使用MEX函数。这种方法的速度要快几个数量级。

以下是代码，供参考：

#include "mex.h"
#include "matrix.h"

/*=======================
 * Output arguments
 *=======================
 */
#define OUT_zCross   plhs[0]
#define OUT_sums     plhs[1]
#define OUT_maxes    plhs[2]

/*=======================
 * Input arguments
 *=======================
 */
#define IN_x        prhs[0]
#define IN_fs       prhs[1]

#define myMax(x,y)     ( ( x ) > ( y ) ? ( x ) : ( y ) )

/*=======================
 * Main Function
 *=======================
 */
void mexFunction ( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
{
    /* params: signal vector; 
       outputs: indices (one-based) of the zero crossings,
       sum of positive values between each pair of zero crossings, and the indices
       (one-based) of the maximum element between each pair
    */

    double *x = NULL;
    double *zCross = NULL;
    double *sums = NULL;
    double *maxes = NULL;
    double curMax = 0;
    unsigned int curMaxPos = 0;
    int Fs = 0;
    unsigned int nZeroCrossings = 0;
    unsigned int nPeaks = 0;
    unsigned int nSamples = 0;
    unsigned int i = 0, j = 0, t = 0;
    bool bIgnoreFirst = false;
    bool bSum = false;

    // Get signal and its size
    x = mxGetPr(IN_x);
    i = mxGetN (IN_x); 
    j = mxGetM (IN_x);

    if (i>1 && j>1) {
        mexPrintf ( "??? Input x must be a vector.\n" );
        return;
    }

    // Length of vector
    nSamples = myMax (i, j);

    zCross = mxCalloc(nSamples, sizeof(double));
    sums = mxCalloc(nSamples, sizeof(double));
    maxes = mxCalloc(nSamples, sizeof(double));

    if (x[0] > 0)
    {
        /* If the first sample is positive, we cannot generate a spike for the first
         pair of zero crossings as this represents part of the signal that is
         negative; therefore, skip the first zero crossing and begin pairing from
         the second */
        bIgnoreFirst = true;
    }
    else if (x[0] == 0)
    {
        // Begin summation from first element
        bSum = true;

        nZeroCrossings = 1;
        sums[0] = x[0];
        curMax = x[0];
        curMaxPos = 0;
    }

    for (t = 1; t < nSamples; ++t)
    {
        // Look for a zero-crossing
        if (x[t] * x[t-1] < 0 || (x[t] == 0 && x[t-1] != 0))
        {
            bool bIgnore = false;

            // If not the first one, we can safely flip the boolean flag
            if (nZeroCrossings != 0)
            {
                bSum = !bSum;
            }
            else if (!bIgnoreFirst)
            {
                // If not, make sure we're not supposed to ignore the first one
                bSum = true;
            }
            else
            {
                bIgnore = true;
            }

            // Store the zero-crossing index
            zCross[nZeroCrossings] = t+1;

            // If this crossing terminated the summation, store and reset the position of the max. element
            if (!bSum && !bIgnore)
            {
                maxes[nPeaks] = curMaxPos+1;
                curMax = 0;
                curMaxPos = 0;
                ++nPeaks;
            }

            ++nZeroCrossings;
        }

        if (bSum)
        {
            sums[nPeaks] += x[t];
            if (x[t] > curMax)
            {
                curMax = x[t];
                curMaxPos = t;
            }
        }
    }

    // Allocate outputs
    OUT_zCross = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);

    OUT_sums = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);

    OUT_maxes = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);

    mxSetPr(OUT_zCross, zCross);
    mxSetM(OUT_zCross, nZeroCrossings);
    mxSetN(OUT_zCross, 1);

    mxSetPr(OUT_sums, sums);
    mxSetM(OUT_sums, nPeaks);
    mxSetN(OUT_sums, 1);

    mxSetPr(OUT_maxes, maxes);
    mxSetM(OUT_maxes, nPeaks);
    mxSetN(OUT_maxes, 1);

    return;
}