Adjusting peaks found in matlab / octave signal - speed issue

Question

I have a function that finds and adjusts peaks based on the rest of the signal. The problem I have is if the signal has 8000 points it works great and is very fast and finishes in about 1min but if it has 200000 points it takes an hour+. Any ideas how to speed up this function?

%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: I'm using octave 3.8.1

I did the profiler like someone recommended but I'm still at a lost as how to improve the speed of this function

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

I think that successive calls to find for logical arrays up tp 200000 elements is a bit waste of resources. Here's an implementation (tested in MATLAB) that should be much faster than the original, even if it uses loops: it passes the array of values just once, no find calls, and adjusts the peaks on the fly:

    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

Please note that %// is just a trick to display correctly the comments on StackOverflow code formatter. otherwise I'd use the normal comment %.

Answer 2

If you sort Af first and save the indexes so that you can later return peak_sig_adj to match the original Af shape, I think you can achieve faster run times as your find() can then be changed to: Faster version of find for sorted vectors (MATLAB)

I haven't done the profiling, but I do think that find() time increases linearly with the size of the input. So you want to use a find function that does not increase so rapidly (e.g. binary search).

Answer 3

If you only ever change the values that were peaks, then you don't need to find on the full matrix except for the first time. There's no way as far as I can see that a part of your signal that didn't originally have a peak value can acquire one.

So we do this once:

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.

The size of peaks therefore decreases with iteration. From a very rough check on much simpler code, this should be faster (how much depends on the ratio of peaks to the actual size of the full signal).