2

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

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

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

当前的实现从实际的数据向量中填充这些列，这样我们之后就不必使用 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 的分配更改为标准 zeros(...) 而不是稀疏数组，那么由于显而易见的原因，缓慢的部分就会变成 sum(...) 和 max(...) 计算。

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

4

1 回答 1

0

现在，我创建了一个 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;
}
于 2012-11-29T09:38:39.283 回答