1-D Bilateral Kernel Filter ... with better than O(N^2) complexity?

I am just testing Bilateral kernel filter for smoothing of 1-D signals. The basic idea of this filtering method is described by the following naive code:

function filteredSignal = bilateralFilter(signal, sigma_spatial, sigma_range) 
% Calculate the size of the 1-D signal 
N = length(signal);
% Initialize the output signal 
filteredSignal = zeros(size(signal));
% Apply bilateral filtering 
for i = 1:N 
    weighted_sum = 0; 
    normalization = 0; 
    
    for j = 1:N 
        % Calculate the spatial difference 
        spatial_diff = abs(i - j); 
        
        % Calculate the range (intensity) difference 
        range_diff = abs(signal(i) - signal(j)); 
        
        % Calculate the bilateral filter weight 
        weight = exp(-spatial_diff^2 / (2 * sigma_spatial^2) - range_diff^2 / (2 * sigma_range^2)); 
        
        % Update the weighted sum and normalization factor 
        weighted_sum = weighted_sum + weight * signal(j); 
        normalization = normalization + weight; 
    end 
    
    % Calculate the filtered value for the current sample 
    filteredSignal(i) = weighted_sum / normalization; 
end 
end 

where

• signal is the input 1D signal you want to filter.
• sigma_spatial and sigma_range are parameters that control the spatial and range domains of the bilateral filter.

The choice of these sigma values depends on the characteristics of your signal and the specific filtering requirements.

My problem is, that for bigger lengths of input signal (N > 1e4) is this naive implementation unacceptable slow (due to the numerical complexity O(N^2)). But in literature is mentioned fact, that this algorithm is possible to modified to the O(c*N) or O(1) complexity, and these "fast" implementations are not available in MATLAB.

I just try to eliminate for-loop over "j", by vectorization, like:

function filteredSignal = bilateralFilter_vec(signal, sigma_spatial, sigma_range) 
% Calculate the size of the 1-D signal 
N = length(signal);
% Initialize the output signal 
filteredSignal = zeros(size(signal));
% Apply bilateral filtering 
for i = 1:N 
    spatial_diff = i-1:-1:i-N; 
    range_diff = abs(signal(i) - signal); 
    weight = exp(-spatial_diff.^2 ./ (2 * sigma_spatial.^2) - range_diff.^2 ./ (2 * sigma_range.^2)); 
    weighted_sum = sum(weight .* signal); 
    normalization = sum(weight); 
    % Calculate the filtered value for the current sample 
    filteredSignal(i) = weighted_sum / normalization; 
end 
end 

The speed-up is significant but still not acceptable.

I will be very happy for any help to transform this code to the more optimized version with better complexity than O(N^2) and still acceptable memory requirements due to the fact, that minimal length of input signal is ~ 1e4 - 1e5 samples.

Answer 1

[Mr. Pavl M.]

Smoothing means LPFKey is in FFT Cooley-Tookey and FFTW (decomposition, recursive devide-and-conqure) realization, did you mean O(N*log(N)) (it can be space-time tradeoff, which is memory storage(RAM) vs computing cycles (clocksteps) in embedded, likely the maximal gain on computation steps saving/minimization, while what is the storage complexity?) complexity? What is this? I work, but no one hired me. Let's amend, employ me.

The main question is what are the feasible high-value, high-utility industrial applications of datum science, control systems engineering and dynamical modes and other specific domain strategies generation to sell the algorithms (find sponsors, buyers, investors, entrepreneours) and their implementations in TCE Matlab, Octave, C, C++, Scheme, Python as the only driving motivation for more in-depth analysis and elaboration within the scope.

This version doesn't run faster, but accurate in TCE GNU Octave (why? because of memory storage enlargement or which mistake to correct / improvement found in it?), while it contains less loops and it smoothes the input:

function filteredSignal = bilFil_vec(signal, sigma_spatial, sigma_range,ver)

%Roll Ver.

% Calculate the size of the 1-D signal
N = length(signal);

% Initialize the output signal
filteredSignal = zeros(size(signal));

% Apply bilateral filtering:

%ver1 (further rapidized as per loop operations minimization):
if ver == 1
    %1. Construct spatial_diff and range_diff matrices as concatenated parametrized growth vectors2d
    %, if not by loop then by reshape or repmat like procedures, how?
    
    %spatial_diff = [0:-1:-N,
    %range_diff = abs(circshift(signal) - signal);
    
    spatial_diff = abs(0:-1:1-N);
    range_diff = abs(signal(1) - signal);
    
    for i = 2:N
        spatial_diff = [spatial_diff;i-1:-1:i-N];
        range_diff = [range_diff;abs(signal(i) - signal)];
    end
    
    %spatial_diff = spatial_diff';
    %range_diff = range_diff';
    
    %replace 2 loops by matrix multiplications:
    weight = exp(-(spatial_diff.^2)./(2*sigma_spatial.^2)-(range_diff.^2)./(2*sigma_range.^2));
    %weighted_sum = sum(weight.*signal);
    %normaliz = abs(weight); %optional
    %normaliz = sum(weight);
    % Calculate the filtered value for the current sample
    %filteredSignal = weighted_sum./normaliz; %*inv(normaliz);
    filteredSignal = signal*weight; %./sum(weight);
else
    %ver 2
    for i = 1:N
        spatial_diff = i-1:-1:i-N;
        range_diff = abs(signal(i)-signal);
        weight = exp(-spatial_diff.^2./(2*sigma_spatial^2)-range_diff.^2./(2*sigma_range^2));
        filteredSignal(i) = sum(weight.*signal)./sum(weight);
    end
end

end

Input signal to it:

Invokation:

filteredSignal = bilFil_vec(inp,2,10,1)

Output produced (as you can see it has been indeed smoothed by my version, which contains internal LPF) (ver == 1):

Please let me know.