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
前述したように、これは遅く、一度に 1 つのチャネルでしか機能しません。遅い部分は (当然のことながら) ループです。行列 G の割り当てをスパース配列ではなく標準の zeros(...) に変更すると、明らかな理由から、遅い部分は sum(...) および max(...) の計算になります。
どうすればこれをより効率的に行うことができますか? 私は MEX 関数を書くことを嫌いではありません。