调整 matlab / octave 信号中的峰值 - 速度问题

Question

我有一个函数可以根据信号的其余部分查找和调整峰值。我遇到的问题是，如果信号有 8000 点，它工作得很好并且非常快并且在大约 1 分钟内完成，但如果它有 200000 点，则需要一个小时以上。任何想法如何加快这个功能？

%example [peaks_found,peak_sig_adj]= rtadjpeak (signal, 3,3)
function [peaks_found,peak_sig_adj] =rtadjpeak (signal, T_multi_lower, T_multi_upper)

    Af=signal;

    % Thresholds
    T_lower = -mean(abs(signal))*T_multi_lower
    T_upper =mean(abs(signal))*T_multi_upper


    % initialisation
    [peaks_r,peaks_c] = find( Af < T_lower | Af > T_upper);
    peaks = find( Af < T_lower | Af > T_upper);
    %find the mean of all the peaks

    counter = 0;

    while ~isempty(peaks)
            peaks = find( Af < T_lower | Af > T_upper);
            try
                    Af(peaks) = ( Af(peaks-1) + Af(peaks+1) ) / 2;
            catch
                    if peaks(1) == 1
                            Af(1) = 0;
                    else
                            Af(end) = 0;
                    end
            end   
            counter=counter+1;
    end
    peaks_found=counter
    peak_sig_adj=Af;
end

PS：我使用的是 octave 3.8.1

我按照别人推荐的方式做了分析器，但我仍然不知道如何提高此功能的速度

profile on;
rtadjpeak(z_sig_combined_L1, 3, 3);
profile off;

>>>T_lower = -0.50551
>>>T_upper =  0.50551
>>>peaks_found =  1013

profshow (profile ("info"));

   #  Function Attr     Time (s)        Calls
---------------------------------------------
  14      find             0.043         1017
   5  binary <             0.023         1017
   1 rtadjpeak             0.023            1
   6  binary >             0.022         1019
  20  binary |             0.018         1015
  23  binary +             0.002         3036
  22  binary -             0.002         1013
  17  binary /             0.001         1013
  21   isempty             0.000         1014
   8  prefix !             0.000         1016
   2       abs             0.000            2
   3      mean             0.000            2
  16       sum             0.000            2
  25   profile             0.000            1
   9     false             0.000            3
   4    nargin             0.000            7
  13      size             0.000            2
   7 isnumeric             0.000            2
  10 binary ==             0.000            4
  11      true             0.000            2

Answer 1

我认为连续调用find 以获取最多 tp 200000 个元素的逻辑数组有点浪费资源。这是一个实现（在 MATLAB 中测试）应该比原来快得多，即使它使用循环：它只传递值数组 once, no find调用，并实时调整峰值：

    function [peaks_found,peak_sig_adj] =rtadjpeak_fast (signal, T_multi_lower, T_multi_upper)

            %// Thresholds
            T_lower = - mean(abs(signal))*T_multi_lower;
            T_upper =   mean(abs(signal))*T_multi_upper;

            %// Initial conditioning for signal
            Af = signal;

            if Af(1) > T_upper
                    Af(1) = T_upper;
            elseif Af(1) < T_lower
                    Af(1) = T_lower;
            end;
            if Af(end) > T_upper
                    Af(end) = T_upper;
            elseif Af(end) < T_lower
                    Af(end) = T_lower;
            end;

            %// Logical array of peaks
            peaks = (Af < T_lower) | (Af > T_upper);

            %// Find continuous peaks, in signal and replace values with average in
            %// the normal range.
            n_max   = numel(peaks);
            in_peak = false;
            counter = 0;
            for k = 1:n_max
                    if ~in_peak && (peaks(k))
                            % Begin of peak
                            n_begin = k;
                            in_peak = true;

                    elseif in_peak && ~(peaks(k))
                            % End of peak, a sample ago
                            n_end   = k-1;
                            in_peak = false;
                            counter = counter + 1;

                            % Calculate average in between
                            n_length = n_end - n_begin + 1;
                            Af_span  = Af(n_end+1) - Af(n_begin-1);
                            Af(n_begin:n_end) = ...
                                  Af(n_begin-1) ...
                                + (Af_span / (n_length + 1)) * (1:n_length);
                    end;
            end;

            %// Set output
            peaks_found  = counter;
            peak_sig_adj = Af;
    end

请注意%// 只是在 StackOverflow 代码格式化程序上正确显示 cmets 的一个技巧。否则我会使用普通评论%。

Answer 2

如果您只更改峰值的值，那么您不需要在完整矩阵上使用find，除非是第一次。据我所知，您的信号中原本没有峰值的部分不可能获得峰值。

所以我们这样做一次：

peaks = find( Af < T_lower | Af > T_upper);

Then inside the loop:

new_peaks = find( Af(peaks) < T_lower | Af(peaks) > T_upper);
peaks = peaks(new_peaks); %making sure indexes match up.

peaks 的大小因此随着迭代而减小。通过对更简单的代码进行非常粗略的检查，这应该会更快（多少取决于peaks 与完整信号的实际大小的比率）。

Answer 3

如果您首先对 Af 进行排序并保存索引，以便稍后返回 peak_sig_adj 以匹配原始 Af 形状，我认为您可以实现更快的运行时间，因为您的 find() 可以更改为： Faster version of find for sorted vectors (MATLAB)

我还没有进行分析，但我确实认为find() 时间随着输入的大小线性增加。所以你想使用一个增长不那么快的查找函数（例如二分查找）。