Matlab：高效地產生聽覺尖峯

Question

我想從MATLAB中的採樣信號中產生一系列「聽覺尖峯」，並且迄今爲止我實施的方法很慢。太慢了。

尖峯產生於第二部分B的http://patrec.cs.tu-dortmund.de/pubs/papers/Plinge2010-RNF.pdf：檢測過零點並對每對之間信號的（平方根壓縮的）正部分進行求和。這給每個釘的高度。他們的位置是通過識別每對過零點之間的最大值來確定樣本的。

我想爲此使用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的分配更改爲標準零（...）而不是稀疏數組，則慢速部分將成爲總和（...）和最大（...）計算，原因很明顯。

我該如何更有效地做到這一點？我不反對寫一個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; 
}