Click here to Skip to main content
14,033,948 members
Click here to Skip to main content
Add your own
alternative version


24 bookmarked
Posted 6 Nov 2014
Licenced CPOL

Impulse noises detection and cancelation in music records

, 6 Nov 2014
Rate this:
Please Sign up or sign in to vote.
Article presents the problem of the detection and cancelation of impulse noise in acoustic signals using an Exponentially Weighted Least Squares (EWLS) and next step error prediction analysis.



When listening to archival music records we often meet (or more precisely hear) the problem of impulse noises occurrence – so called ‘crackles’. It is getting worse when listening to vinyl records or cassette tapes.

There are many algorithms whose task is to remove such distortions – some are suitable for hardware implementation, some of them can be used only in offline implementations.

This article is to present adoption of identification method Exponentially Weighted Least Squares (EWLS) and next step error prediction analysis.

Despite the fact that this algorithm will be used to detect and cancel noise in recorded music, it can successfully be used to cancel impulse noise from any set of even multidimensional signals.

Acoustic signals

Acoustic signals which are present in our real world are time-continuous signals, and in this form are recorded on analog media such as vinyl records or cassettes.

Recording on the more modern digital data media (CDs, hard drives, Flash memory) requires digitization of the analog signal - sampling and quantization.

Sound frequencies heard by the human ear are in the range of about 20 kHz. According to Nyquist - Shannon theorem the sampling frequency should be at least twice of highest frequency peak in signal spectrum in order to avoid the effect of aliasing. In the processing of music we used set sampling frequency of 44.1 kHz.

In the article it will be assumed that processed signal is digital and was sampled correctly according to Nyquist – Shannon theorem.

Signal estimation

Digital acoustic signal can be consider as discrete time function f(t), where t = 0, 1, … represents normalized time. Let assume that next function values (signal samples) can be represented by autoregressive model (AR):

where  represents signal sample value in moment ,  represents time shift operator, represents rank of AR model, represents ith model parameter,  represents white noise.

Let us assume vectors:


Using them we are able to reorder autoregressive model equation:

If we will assume that vector  represents real (unknown to us) parameters of  signal, then we can assume that there is a vector  which is estimation of vector .

Recursive identification algorithms calculate prediction of next signal value  according to data (values) collected till this step:


Recursive identification algorithm

Due to equation above, if want to calculate signal value prediction for the next step we need to designate vector , as vector  is given as an outcome of previous operations.

To calculate  we will use Exponentially Weighed Least Squares algorithm. EWLS algorithm, as every in the least squares algorithms family, calculates model parameters by minimizing cost function on least squares – in EWLS it gives exponentially less weight to previous error values.

To reduce the impact of older values to the current calculation we need to assume a so called forgetting factor . The forgetting factor takes value from 0 to 1, wherein if  then EWLS algorithm’s memory is shorter.

The recursive EWLS is given with equation:

where  is an error function – the difference between true signal value in  moment and estimated signal value in  moment:

and  is a gain vector , described as:

where matrix  can be calculated with equation:


Interference detection

Impulse noises in two dimensional signal can be seen as samples with amplitude much higher than neighboring samples (Fig.1).

Fig.1. Impulse noise example in acoustic signal

To detect such a samples we will use algorithm based error function prediction, defined as

Let us declare two additional parameters, the first one is prediction error local variance:

the second one is standard deviation of prediction error:

Noise detection refers to verify the value of . Naturally, higher value of error function in  moment indicates that signal value in this moment is disturbed.

Obviously, the condition given above is strongly fuzzy – what does it mean ‘higher value of error function’?

To complete this condition we assume that suspicious samples of acoustic signal induce error function values above scaled standard deviation of this function:

where  is scaling coefficient which is responsible for algorithm sensitivity and need to be defined as algorithm input data. Note that too high  value may cause missing disturbed samples, whereas too low value may cause classifying ‘clean’ samples as disturbed ones.

Denote by  binary function with possible values 0 (if tth sample is ‘clean’) and 1 (if tth value is disturbed). Then:

Error prediction local variance can be calculated using:

with initial condition .


Reconstruction of distrubed samples

Sample marked as disturbed is to  be deleted and replaced with reconstructed one. For this purpose we will apply linear interpolation using neighboring ‘clean’ samples.

Another issue to be consider is the maximum number of next samples that can be disturbed. Due to this different  parameter value should be places in equation below

where indicates position of sample under reconstruction in the step.



Initial parameter were taken in algorithm simulation:

A 20 second length music record was processed with algorithm, results can be seen at Fig. 2 and 3.

Fig.2. Simulation results. Upper plot presents original (disturbed) signal, lower plot presents this signal after detecting and canceling noises

Fig.3. Zoom of dataset from Fig. 2


Use of code

Implementation of this algorithm in Matlab environment is quite easy, as a lot of predefined functions can support us.

Use of code is by calling function:

impulse_noise_elimination(file_1, file_2, model_order, n, lambda, p, save_ar)


%Function Impulse Noise Elimination input parametres:
% * file_1 is a name of input track, ex. 'music.wav';
% * file_2 is a name of output track, ex. 'music_mod.wav';
% * model_order is a choosen number of estimated model parameters; can be
%   choosen from range: 1 <= model_order <= 25; mostly the best choice is
%   model_order = 4;
% * n is a scaling coefficient of local variance in detecting algorithm;
%   its value decides of the level of impulse noise samples; mostly n = 3
%   is the best choice, but n is from range 2 <= n <= 20;
% * lambda is a memory coefficient in EW - LS algorithm; its value
%   dependends on track dynamic; lambda is from the range: 
%   0.001 (low dynamic) <= lambda <= 0.999 (high dynamic);
% * p is a scaling coefficient of P matrix (eye matrix rank model_order)
%   starting value; p is from range 10^3 <= p <= 10^6;
% * save_ar = 0 or 1; 0 if model parametres should not be saved, 1 if model
%   parameters should be saved.

Milestones of this function are shown below.

Input data check

[track fs tmp] = wavread(file_1); %orginal track (wav file)


if model_order > 25
    model_order = 25;
elseif model_order < 1
    model_order = 1;

if n > 20
    n = 20;
elseif n < 2
    n = 1;

if lambda > 1
    lambda = 0.999;
elseif lambda < 0
    lambda = 0.001;

if p < 1000
    p = 1000;
elseif p > 1000000
    p = 1000000;

if save_ar ~= 1
    save_ar = 0;
elseif save_ar == 1
    model_vector = zeros(max(size(track)), model_order);

Main loop steps

Initial parameters:

%Starting parametres of EW - LS algorithm:
org_track = track;
model = zeros(model_order, 1);
K = 0;
P = p * eye(model_order);
variance = 0;

Loop definition:

for t = (model_order + 1) : (max(size(track)) - 4)

% All operations


As it can be seen, loop start condition is

(model_order + 1)

wheras the end condition is 

(max(size(track)) - 4)

as we need to collect at least model rank samples number to start identification and our implemention considers that next four samples can be distrubed.

To get another signal value and calculate prediction error we do:

signal = track(t-1 : -1 : t-model_order);
prediction_error = track(t) - model' * signal;

When previous sample was 'clean', variance can be calculated:

variance = lambda*variance + (1 - lambda)*(prediction_error)^2;

Otherwise variance calculation does not change. To initiate we use:

variance = prediction_error^2;

The simplest (non-optimal) way to check next four samples (when first one is disturbed) is by using if condition:

if abs(prediction_error) > n * sqrt(variance) %detecting test
    noise_samples = 1;
    signal_1 = track(t : -1 : t-model_order+1);
    prediction_error_1 = track(t+1) - model' * signal_1;

    if abs(prediction_error_1) > n * sqrt(variance)
        noise_samples = 2;
        signal_2 = track(t+1 : -1 : t-model_order+2);
        prediction_error_2 = track(t+2) - model' * signal_2;
        if abs(prediction_error_2) > n * sqrt(variance)
            noise_samples = 3;
            signal_3 = track(t+2 : -1 : t-model_order+3);
            prediction_error_3 = track(t+3) - model' * signal_3;
            if abs(prediction_error_3) > n * sqrt(variance)
                noise_samples = 4;

    noise_samples = 0;
    variance = lambda*variance + (1 - lambda)*(prediction_error)^2;

After that we can try to reconstruct disturbed samples:

if noise_samples ~= 0
    for i = 0 : noise_samples - 1
        track(t + i) = track(t-1) + (-1 / (noise_samples + 1))*(track(t-1) - track(t + noise_samples))*(i+1);
    noise_samples = 0;

At the end we can update our estimator parameters

prediction_error = track(t) - model' * signal;
K = (P * signal) / (lambda + signal'*P*signal);
P = (1 / lambda)*(P - (P*signal*signal'*P)/(lambda + signal'*P*signal));
model = model + K * prediction_error;

and start all over again.


Presented recursive algorithm can be implemented in real-time applications. The delay produced by this algorithm is defined by the maximum number of next samples that can be disturbed.

Acoustic signals considered in this article by the nature are characterized by their high dynamic of amplitude changes. Due to this are signals types can successfully be processed with this algorithm (for example measurement data).


[1] Niedźwiecki M. Kłaput T. Fast recursive basis function estimators for identification of time-varying processes., IEEE Transactions on signal processing, vol. 50, no. 8, pp. 1925 - 1934. August, 2002.

[2] Söderström T. Stoica P. System identyfication. Prentice-Hall International, Hemel Hempstead, 1989.

[3] Zieliński Tomasz P. Cyfrowe przetwarzanie sygnałów. Wydawnictwa Komunikacji i Łączności, Warszawa, 2005.




This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


About the Author

Jakub Szymanowski
Systems Engineer Enspirion LLC (ENERGA Group)
Poland Poland
After graduating with a degree of Master of Science, Engineer in Automation and Robotics (at the Gdansk University of Technology) began his career as a Service Engineer of Ships' Systems. He spent four years in this profession and had a lot of places all over the world visited.

In 2012 he started to work for Polish energy company ENERGA (ENERGA Agregator LLC, ENERGA Innowacje LLC and Enspirion LLC) as Systems Engineer and Solutions Architect. Thereafter he participated in few large-scale research projects in the field of Information Technology and Systems Integration for Energy & Power Management and Demand Response applications as a Solution Architect, Business Architect or Project Manager.

His private interests include Signal and Image Processing, Adaptive Algorithms and Theory of Probability.

You may also be interested in...


Comments and Discussions

QuestionVery interesting Pin
Kenneth Haugland16-Apr-15 6:37
professionalKenneth Haugland16-Apr-15 6:37 
QuestionImpulse noise detection and cancelation Pin
Suresh Doraiswamy6-Nov-14 22:17
professionalSuresh Doraiswamy6-Nov-14 22:17 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Cookies | Terms of Use | Mobile
Web01 | 2.8.190423.1 | Last Updated 6 Nov 2014
Article Copyright 2014 by Jakub Szymanowski
Everything else Copyright © CodeProject, 1999-2019
Layout: fixed | fluid