Recursive Norm Computation

Recursive Computation of Signal Vector Norms

written by Gerhard Schmidt, Tim Owe Wisch, and Katharina Rebbe

Problem

In several signal processing applications signal vectors that contain the last \(N\) sample are utilized. Those vectors are usually defined as

\begin{equation} {\bf x}(n) \,\,=\,\, \Big[ x(n),\,x(n-1),\,x(n-2),\, ...,\,x(n-N+1)\Big]^{\textrm{T}}. \label{eq_red_rec_comp_vec_norms_00} \end{equation}

Furthermore, some applications require to compute the squared norm of such vectors. A direct computation according to

\begin{equation} \big\| {\bf x}(n)\big\|^2 \,\,=\,\, \sum\limits_{i=0}^{N-1}x^2(n-i) \label{eq_red_rec_comp_vec_norms_01} \end{equation}

is numerically quite robust, but requires also \(N\) multiplications and \(N-1\) additions every sample. In order to show the numerical robustness, we generated a signal that contains white Gaussian noise. Every 1000 samples we varied the power by 20 dB (up and down in an alternating fashing). From that signal we extracted signal vectors of length \(N=128\) and computed the squared norm accordign to Eq. \(( \ref{eq_red_rec_comp_vec_norms_01} )\) in floating point precision, once with 64 bits and once with 32 bits and depict the results in a logarithmic fashion in Fig. 1. Additionally, the difference between the two versions is shown in the lowest diagram.



	Fig. 1: Input signal and iterative norm computations.

Recursive Computation

Remark:

The following ideas (and the solutions) can also be used for other recursive computations such as mean estimations \(y(n) = \frac{1}{N}\sum\limits_{i=0}^{N-1}x(n-i)\).

If the norm of the signal vector has to be computed every sample, often recursive computations are favoured since this lead to a significant recution of computational complexity - especially for signal vectors with a large amount of elements [HS 2004]. The recurive variant starts with an initialization. Is is assumed that the signal vector contains zeros at time index \(n=0\) and thus also the squared norm is initialized with zero:

\begin{eqnarray} {\bf x}(0) \,\,&=&\,\, \big[ 0,\,0,\,0,\, ...,\,0\big]^{\textrm{T}}, \\[2mm] \big\| {\bf x}(0)\big\|^2 \,\,&=&\,\, 0. \label{eq_red_rec_comp_vec_norms_02} \end{eqnarray}

Since with every sample only one new sample value is added to the signal vector and one sample value (the oldest) is leaving the vector, the norm can be computed recursively according to

\begin{equation} \big\| {\bf x}(n)\big\|^2 \,\,=\,\, \big\| {\bf x}(n-1)\big\|^2 + x^2(n) - x^2(n-N).\label{eq_red_rec_comp_vec_norms_03} \end{equation}

Using this "trick" only two multiplications and two additions are required to update the norm. However, from a numerical point-of-view, this computation is not as robust as the direct approach according to Eq. \(( \ref{eq_red_rec_comp_vec_norms_01} )\).



	Fig. 2: Input signal and recursive as well as iterative norm computations.

The problem with the recursive computation according to Eq. \(( \ref{eq_red_rec_comp_vec_norms_03} )\) is that a squared signal value \(x^2(n)\) is used twice in the update rule:

once when it enters the signal vector and
once when it leaves the vector.

If the computation is done in floating-point arithmetic first the mantissa of the values that should be added or subtracted are adjusted (shifted) such that the exponents of both values are equal. This can be interpreted as some sort of quantization. If the norm value has changed between the "entering" and the "leaving" event, a small error occurs. Unfortunately, this error is biased and thus error accumulation appears.

Usually, this is not critical, because the error is really small, but if the signal has large power variations and the recursion is performed several thousand times, the small error might become rather large.

Due to that problem the norm might get negative, which leads - in some cases - to severe probelms. E.g. several gradient based optimizations perform a division by the norm of the excitation vector. If the norm gets negative it means that the direction of the gradient is switched and divergence might be the result. This could be avoided by limiting the result of the recursive computation by the value 0:

\begin{equation} \big\| {\bf x}(n)\big\|^2 \,\,=\,\, \max\Big\{0,\,\big\| {\bf x}(n-1)\big\|^2 + x^2(n) - x^2(n-N)\Big\}. \label{eq_red_rec_comp_vec_norms_04} \end{equation}

This improves robustness, but does not help against error accumulation as depicted in Fig. 2.

Mixed Recursive/Iterative Computation

A solution to this error accumulation problem is the extent the recursive computation according to Eq. \(( \ref{eq_red_rec_comp_vec_norms_04} )\) by an iterative approach that "refreshes" the recursive update from time to time. This can be realized by adding in a separate variable \(N_{\textrm{rec}}(n)\) all squared input samples:

\begin{equation} N_{\textrm{rec}}(n) \,\,=\,\, \left\{ \begin{array}{ll} x^2(n), & \textrm{if}\,\mod(n,N)\,\equiv\, 0, \\[2mm] N_{\textrm{rec}}(n-1) + x^2(n), & \textrm{else}. \end{array}\right. \label{eq_red_rec_comp_vec_norms_05} \end{equation}

If \(N\) samples are added this variable is replacing to the recursively computed norm and the original sum \(N_{\textrm{rec}}(n)\) is reinitialized with 0:

\begin{equation} \big\| {\bf x}(n)\big\|^2 \,\,=\,\, \left\{ \begin{array}{ll} N_{\textrm{rec}}(n), & \textrm{if}\,\mod(n,N)\,\equiv\, N-1, \\[2mm] \max\Big\{0,\,\big\| {\bf x}(n-1)\big\|^2 + x^2(n) - x^2(n-N)\Big\},& \textrm{else}. \end{array}\right. \label{eq_red_rec_comp_vec_norms_06} \end{equation}

The additional mechanism adds only a few additions, but helps a lot against error accumulation as indicated in the last example of this section depicted in Fig. 3. Thus, if you face problems with norms of signal vectors you might think about using this mixed method.



	Fig. 3: Input signal and mixed recursive/iterative as well as purely iterative norm computations.

References

[HS 2004]

E. Hänsler, G. Schmidt: Acoustic Echo and Noise Control, Wiley, 2004.

Code Example

As usual you can find below some matlab code that shows examples on the inidvidual methods of norm computations as well as the corresponding analyses.

{slider title="Different methods of computing vector norms" class="grey" open="false"}

%**************************************************************************
% Parameters
%**************************************************************************
N            =   128;
Sig_duration = 10000;

%**************************************************************************
% Generate input signal
%**************************************************************************

% Generate white Gaussian noise
sig = single(randn(Sig_duration,1));

% Boost every second 1000 signal values by 60 dB  
for k = 1001:2000:Sig_duration;
    sig(k-1000:k) = sig(k-1000:k) * 1000;
end;

%**************************************************************************
% Compute norms
%**************************************************************************
x_vec               = single(zeros(N,1));
Norm_rec_curr       = single(0);
Norm_rec_curr_lim   = single(0);
N_rec               = single(0);
C_rec               = single(0);
Norm_mixed_curr     = single(0);
N_rec_lim           = single(0);
C_rec_lim           = single(0);
Norm_mixed_lim_curr = single(0);

Norm_double_prec    = zeros(Sig_duration,1);
Norm_iterative      = single(zeros(Sig_duration,1));
Norm_recursive      = single(zeros(Sig_duration,1));
Norm_recursive_lim  = single(zeros(Sig_duration,1));
Norm_mixed          = single(zeros(Sig_duration,1));
Norm_mixed_lim      = single(zeros(Sig_duration,1));

for k = 1:Sig_duration

    %**********************************************************************
    % Update signal vector (not very efficient, but o.k. for here
    %**********************************************************************
    x_new      = sig(k);
    x_old      = x_vec(N);
    x_vec(2:N) = x_vec(1:N-1);
    x_vec(1)   = x_new;
    
    %**********************************************************************
    % Norm in double precicion
    %**********************************************************************
    Norm_double_prec(k) = double(x_vec)' * double(x_vec);
    
    %**********************************************************************
    % Iterative norm (of course, there exixt optimized Matlab functions
    % for this purpose, but that's another "story")
    %**********************************************************************
    for n = 1:N
        Norm_iterative(k) = Norm_iterative(k) + x_vec(n)*x_vec(n);
    end;
    
    %**********************************************************************
    % Recursive norm computation
    %**********************************************************************
    Norm_rec_curr     = Norm_rec_curr + x_new^2 - x_old^2;
    Norm_recursive(k) = Norm_rec_curr;
    
    %**********************************************************************
    % Recursive norm computation
    %**********************************************************************
    Norm_rec_curr_lim     = Norm_rec_curr_lim + x_new^2 - x_old^2;
    Norm_rec_curr_lim     = max(0, Norm_rec_curr_lim);
    Norm_recursive_lim(k) = Norm_rec_curr_lim;  
    
    %**********************************************************************
    % Mixed compuation of the norm
    %**********************************************************************
    Norm_mixed_curr = Norm_mixed_curr + x_new^2 - x_old^2;
    
    C_rec = C_rec + 1;    
    if (C_rec == N)
        C_rec = 0;
    end;
   
    if (C_rec == 0)
        N_rec = 0;
    end;
    N_rec = N_rec + + x_new^2;
    
    if (C_rec == N-1)
        Norm_mixed_curr = N_rec;    
    end;
    
    Norm_mixed(k) = Norm_mixed_curr; 
    
    %**********************************************************************
    % Mixed compuation of the norm with limiation
    %**********************************************************************
    Norm_mixed_lim_curr = Norm_mixed_lim_curr + x_new^2 - x_old^2;
    Norm_mixed_lim_curr = max(0,Norm_mixed_lim_curr);
    
    C_rec_lim = C_rec_lim + 1;    
    if (C_rec_lim == N)
        C_rec_lim = 0;
    end;
   
    if (C_rec_lim == 0)
        N_rec_lim = 0;
    end;
    N_rec_lim = N_rec_lim + + x_new^2;
    
    if (C_rec_lim == N-1)
        Norm_mixed_lim_curr = N_rec_lim;    
    end;
    
    Norm_mixed_lim(k) = Norm_mixed_lim_curr;
    
end;

%**************************************************************************
% Show results
%**************************************************************************
fig = figure(1);
set(fig,'Units','Normalized');
set(fig,'Position',[0.1 0.1 0.8 0.8]);

t = 0:Sig_duration-1;

subplot('Position',[0.07 0.8 0.9 0.17]);
plot(t,10*log10(Norm_double_prec),'b', ...
     t,10*log10(max(0.01, Norm_iterative))+1,'r');
grid on
set(gca,'XTickLabel','');
ylabel('dB')
legend('Norm in double precision', ...
       'Iteratively computed norm (lifted by 1 dB)');
axis([0 Sig_duration-1 0 120]);

subplot('Position',[0.07 0.62 0.9 0.17]);
plot(t,10*log10(Norm_double_prec),'b', ...
     t,10*log10(max(0.01, Norm_recursive))+1,'r');
grid on
set(gca,'XTickLabel','');
ylabel('dB')
legend('Norm in double precision', ...
       'Recursively computed computed norm (lifted by 1 dB)');  
axis([0 Sig_duration-1 0 120]);   

subplot('Position',[0.07 0.44 0.9 0.17]);
plot(t,10*log10(Norm_double_prec),'b', ...
     t,10*log10(max(0.01, Norm_recursive_lim))+1,'r');
grid on
set(gca,'XTickLabel','');
ylabel('dB')
legend('Norm in double precision', ...
       'Recursively computed computed norm with limitation (lifted by 1 dB)');     
axis([0 Sig_duration-1 0 120]);   
   
subplot('Position',[0.07 0.26 0.9 0.17]);
plot(t,10*log10(Norm_double_prec),'b', ...
     t,10*log10(max(0.01, Norm_mixed))+1,'r');
grid on
set(gca,'XTickLabel','');
ylabel('dB')
legend('Norm in double precision', ...
       'Mixed recursively/iteratively computed norm (lifted by 1 dB)');     
axis([0 Sig_duration-1 0 120]);   
   
subplot('Position',[0.07 0.08 0.9 0.17]);
plot(t,10*log10(Norm_double_prec),'b', ...
     t,10*log10(max(0.01, Norm_mixed_lim))+1,'r');
grid on
xlabel('Samples');
ylabel('dB')
legend('Norm in double precision', ...
       'Mixed recursively/iteratively computed norm with limitation (lifted by 1 dB)');   
axis([0 Sig_duration-1 0 120]);

{/sliders}

Basic Idea of RED

RED (robust and efficient DSP) is a collection of algorithms that should help engineers to implement basic signal processign tasks in - as the name states - a robust and efficient way. This includes mainly numerical aspects. An overview about the RED collection can be found here.

Here you can find the chapter on robust recursive norm computations as a single pdf file. If you would like to download all chapters please go to the RED overview page and get the pdf document that contain all chapters from there.

Authors of this Section:

Katharina Rebbe
Gerhard Schmidt
Owe Wisch

Related Downloads:

	Pdf file of this section

	Matlab script for a computing vector norms in an efficient and robust manner.

We are delighted to announce that Marten Finck has been awarded the GaS prize for his Master's thesis on "Evaluation of Speech Quality Using Machine Learning".

In his work, he investigated possible side effects of magnetic resonance-guided focussed ultrasound treatment (MRgFUS) on speech quality in patients with neurological diseases such as Parkinson's disease and essential tremor. This innovative and non-invasive treatment method is used to reduce tremor, i.e. the shaking of limbs.

To evaluate the effects on speech, pre- and post-operative speech recordings of 56 patients were analysed: Speech quality was assessed in five categories (monotony, clarity, fluency, auditory effort and vocal tremor) through a survey with 21 participants and used to create a convolutional neural network (CNN). This generates an objective assessment of speech quality using mel frequency cepstral coefficients (MFCCs), fundamental frequency and loudness. A mean absolute error (MAE) of 0.81 on a scale of 0 to 10 was achieved, compared to a standard deviation of 0.67 for the survey participants compared to a reference assessment. The evaluation of the survey revealed minor positive and negative changes in speech quality, which were balanced overall and tended to indicate natural speech variations in connection with the patients' underlying illnesses.

Congratulations on this well-deserved honour!