Audioread Impulse noise detection and removal of impulse signals

Please find the attached source code file and audio file * .wav file.

I would like to filter out the impulse noise from the right channel.

Answer 1

[Adam Danz]

This uses findpeaks to locate the impuse signal in channel 2 of the wave file. Peaks are required to be at least 3 standard deviations larger than the mean signal. The values are replaced with the median of that channel (excluding the impulse signal). fFilt is the filtered version of f.

% Read in wav file
FileName = 'rec_noise.wav';  %full path is better
[f,fs]    = audioread(FileName);

% Locate peaks
minHeight = mean(abs(f(:,2))) + std(abs(f(:,2)))*3;  % min height to be 'peak';  3 std above mean
[pks, locs] = findpeaks(abs(f(:,2)), 'MinPeakHeight',minHeight);
pks = pks .* sign(f(locs,2)); % pks is only used for plotting

% Replace peaks with the median value of f(:,2) expluding outliers.
% You could also just remove those values but you'd need to remove them in chan 1, too.
isOut = ismember(1:size(f,1), locs); 
replacementValue = median(f(~isOut,2));  % you could also use NaN, 0, etc....
fFilt = f; 
fFilt(isOut,2) = replacementValue; 

Plot the results.

clf()
% Plot channel 1
sp(1) = subplot(3,1,1);
plot(sp(1), f(:,1),'DisplayName','chan 1')
title(sp(1), 'channel 1')

% Plot channel 2
sp(2) = subplot(3,1,2);
plot(sp(2), f(:,2),'DisplayName','chan 2')
title(sp(2), 'channel 2')
% Show outliers (impulses)
hold(sp(2),'on')
plot(sp(2),locs, pks, 'rx', 'DisplayName', 'Outliers')

% Plot the filtered channel 2
sp(3) = subplot(3,1,3);
plot(sp(3), fFilt(:,2),'DisplayName','chan 2 Filtered')
title(sp(3), 'channel 2 filtered')
linkaxes(sp)

[Jogger]

Thanks a lot .. super solution.

[Jogger]

I tested the solution which works well with all impulse noise . But I don't know sometimes - I am getting following error

Error using findpeaks
Too many input arguments.

Error in AudioRecording (line 50)
[pks, locs] = findpeaks(abs(f(:,2)));

[Adam Danz]

Matlab, unrelated to your error, I made a change to my answer because the older version was identifying peaks that weren't really peaks and therefore removing good data.

The updated version now requires that a peak be at least 3 standard deviations from the mean signal. As you can see, there are much fewerd red x markers in the 2nd subplot.

The changes were

% Old
[pks, locs] = findpeaks(abs(f(:,2))); 

% New
minHeight = mean(abs(f(:,2))) + std(abs(f(:,2)))*3;  % min height to be 'peak';  3 std above mean
[pks, locs] = findpeaks(abs(f(:,2)), 'MinPeakHeight',minHeight);

Regaring the error you reported, the older version only had 1 input and the findpeaks function allows for up to two inputs plus name-value property pairs. So something's fishy.

1. Check that you don't have a different function named "findpeaks". After you experience the error, run this line of code and share the results: which('findpeaks', '-all')
2. Are you using the exact same variable f from the AudioRecording.m file?

[Jogger]

• Check that you don't have a different function named "findpeaks". After you experience the error, run this line of code and share the results: which('findpeaks', '-all')

Please find the output

K>> which('findpeaks', '-all')
C:\Users\shastry\Release_version\Audio\findpeaks.m
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\@dspdata\findpeaks.m  % Shadowed dspdata method
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\findpeaks.m           % Shadowed

• Are you using the exact same variable f from the AudioRecording.m file?

Yes, I am using the same ..

 [pks, locs] = findpeaks(abs(f(:,2)), 'MinPeakHeight',minHeight);
K>> minHeight = mean(abs(f(:,2))) + std(abs(f(:,2)))*3;  % min height to be 'peak';  3 std above mean
K>> minHeight

minHeight =

    0.0325

K>> [pks, locs] = findpeaks(abs(f(:,2)), 'MinPeakHeight',minHeight);
Error using findpeaks
Too many input arguments.

K>> which('findpeaks', '-all')
C:\Users\shastry\Release_version\Audio\findpeaks.m
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\@dspdata\findpeaks.m  % Shadowed dspdata method
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\findpeaks.m           % Shadowed

In addition - would like here some more data plotting like, Signal Amplitude, Frequency and phase information.

updated version now requires that a peak be at least 3 standard deviations from the mean signal

It would be helpful for me if you can provide a small concept , as I have few more *.wav file to analyze.

Thank you!

[Adam Danz]

There's the problem. You're not using Matlab's findpeaks function. The "findpeaks" function listed at the top must be a custom function and since it has the same name and it's on path, it's shadowing Matlab's function.

I suggest you rename the file and the function (if it's a function). Then run the "which" command again to make sure the problem is fixed. For example, this is what I see when I run it.

>> which findpeaks -all
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\findpeaks.m
C:\Program Files\MATLAB\R2019b\toolbox\signal\signal\@dspdata\findpeaks.m  % dspdata method

I'm not sure what you're asking about regarding the 'small concept'. The findpeaks function can be very sensitive and sometimes it identifies peaks that aren't really peaks. Prior to making the 3*std() change, too many "peaks" were being identified and you were losing some of the signal. To make findpeaks less sensitive, the new change requires that peaks are at least 3 standard deviations from the mean of the signal.

[Jogger]

function listed at the top must be a custom function!!

Super catch.. My initial preparation was to find the peaks , so I created a function findpeak ...I didn't reliaze toolbox has similar one.

Again thank you very much.

[Jogger]

What if we have to apply filtering for both channels ?

[Jogger]

Please find the attached *.wav file (different from earlier attachment) and my latest test script .

By analyzing the result from the attached file, I think impulse noise is has a mis detection and so not removed.

Can you please check and help me ?

Thank you!

[Adam Danz]

"What if we have to apply filtering for both channels ?"

This line below applies the filtering to channel 2. You can repeat the process for other channels.

[pks, locs] = findpeaks(abs(f(:,2)), 'MinPeakHeight',minHeight);

[Adam Danz]

"I think impulse noise is has a mis detection "

Based on the image you shared, the impulses are detected pretty well. Look at the upper two subplots. But it looks like you're not indexing correctly and therefore, the impulses aren't being removed correctly.

There are just a few lines in my code. Take some time to understand what each one is doing.

This line below is removing the detected impuses from chan 2. Make sure you're accessing the correct channel index number.

fFilt(isOut,2) = replacementValue; 

[Jogger]

Hi Adam,

Please find the code snipet I am using for chan-1.

I see no noise filtering is applied on channel-1 because there is no replace value selected. I saw closely on minHeight1 value . It does not fit accuratley for above attached noise file ( Fail_4.zip).

section below

minHeight1 = mean(abs(f(:,1))) + std(abs(f(:,1)))*3; % min height to be 'peak';

replacementValue1 ==> has a value 0. So no filtering is appiled.

I need help in selecting best algorithm in determining the constant[minHeight1 [-1 fraction +1]] which is best fit all/most impulse noise cases.

[path]      = get_currentdir();
files       = dir(path);

for fileIndex=1:length(files)
    if (files(fileIndex).isdir == 0)
        if (~isempty(strfind(files(fileIndex).name,'wav')))
            FileName           = fullfile(path,files(fileIndex).name);
            [f,fs]             = audioread(FileName);
        end
    end
end
%%
pOrig = audioplayer(f,fs);
pOrig.play
%%
figNo = 200;
clf();figure(figNo+2)

minHeight1 = mean(abs(f(:,1))) + std(abs(f(:,1)))*3;  % min height to be 'peak';  3 std above mean


%-----------------------------------------------------------------------------------%
[pks1, locs1] = findpeaks(abs(f(:,1)), 'MinPeakHeight',minHeight1); % MinPeakHeight,threshold is more suitable for File_4
pks1 = pks1 .* sign(f(locs1,1)); % pks is only used for plotting

%% Plot channel 1
sp(1) = subplot(2,1,1);
plot(sp(1), f(:,1),'DisplayName','chan 1')
title(sp(1), 'channel 1')
% Show outliers (impulses)
hold(sp(1),'on')
plot(sp(1),locs1, pks1, 'g-x', 'DisplayName', 'Outliers')

isOut1 = ismember(1:size(f,1), locs1); 
replacementValue1 = median(f(~isOut1,1));  % you could also use NaN, 0, etc....
title(sprintf('MinPeakHeight RepVal:%4.3f',replacementValue1));

Thanks a lot!

[Jogger]

I have raise a new link for above request

https://www.mathworks.com/matlabcentral/answers/505525-findpeaks-algorithm-to-determine-constant-value-for-impulse-noises

Thank you!

[Adam Danz]

You're not applying the filtering.

You're missing the last two lines from the first block of code in my answer.

The last two lines of your code above are doing the following:

% this next line produces a logica index identifying the outliers/peaks
isOut1 = ismember(1:size(f,1), locs1); 

% this next line merely calculates a value that will be used 
% to replace the peak values.  For example, if the f(100) has a value of 10
% and is identified as a peak, its value will be replaced with replacementValue
% f(100) = replacementValue.  
replacementValue1 = median(f(~isOut1,1));  % you could also use NaN, 0, etc....

"replacementValue1 ==> has a value 0. So no filtering is appiled."

As explained above, a value of 0 is reasonable for this variable and this variable is not applying the filtering.

[Jogger]

You mean below lines were missing ?

fFilt1 = f; 
fFilt1(isOut1,1) = replacementValue1; 

I have applied be lines as well ,

I play back the audio using.

pL = audioplayer(fFilt(isOut1,1), fs);

pL.play;

[Adam Danz]

Yes.

% This line copies the audio data so you don't lose it when you replace the peaks. 
fFilt1 = f; 

% This line does the filtering.  It replaces the peaks (isOUt) with the value stored in replacementValue
fFilt1(isOut1,1) = replacementValue1; 

[Adam Danz]

This line is doing the exact opposite of what the filter is doing.

% In this line you're accepting **only the peaks**
pL = audioplayer(fFilt(isOut1,1), fs);

% In this line, the peaks are being filtered-out
fFilt1(isOut1,1) = replacementValue1; 

This line would be correct, with the "~", but keep in mind, this method reduces the number of samples in f. It completely removes the peak values. In my answer, the peak values are replaced with the mean signal. Either approach is fine - just know the difference.

pL = audioplayer(fFilt(~isOut1,1), fs);

Please take a few minutes to go through each line of my answer and understand exactly what it's doing, what the variables are, their shape and size, their content, etc.

After reading in the file, there's only 7 lines of code in my answer. The rest is just plotting.

[Jogger]

I missed ~ factor . My mistake . Understood . Thanks you very much.

[Adam Danz]

Glad I could help!

[Jogger]

Hi Adam Danz , Guillaume

I have raised a new issue. Please find the below link for more information

https://www.mathworks.com/matlabcentral/answers/507386-remove-noise-in-video

Can you please help me ? Thank you!

[Jogger]

Hi Adam Danz , Guillaume

Can you please help me with below issue in the link ?

https://www.mathworks.com/matlabcentral/answers/511453-ismember-is-not-supported-for-fixed-point-conversion

Thanks !