9
$\begingroup$

I have some underwater recordings (wav files) that have periods of relatively high self noise (something bumping the hydrophone). These have the potential of biasing the ambient sound measures I want to do by increasing the measured sound (RMS de re 1uPa).

I'd like to run a detector which identifies these noise events and then removes them so they are not carried forward when I measure RMS of background sound. The noises are broadband sounds whose bandwidth overlaps with my signals of interest so a band pass filter won't work. I am familiar with matlab if that helps.

Any suggestions?

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Welcome @duck. Could you provide a short spectrogram of the sound+ noise to help get targeted answers. Also what is meant by 'background sound' vs the 'noise'. Typically background and noise are used interchangeably. $\endgroup$
    – Thejasvi
    Commented Sep 8, 2022 at 5:45
  • 1
    $\begingroup$ @Thejasvi I added a picture and some additional context $\endgroup$
    – ducks
    Commented Sep 8, 2022 at 12:58

2 Answers 2

4
$\begingroup$

There is no entirely fail safe method for this, but as I understand your question you want to remove an intermittent short noise (the bumping noise) from an otherwise relatively constant noise level recording.

I have a few suggestions that can be run easily in MatLab/Octave/Python provided your recording is not very long:

As @WMXZ suggests in another response a simple threshold filter might work for you: Find a level which all your "ambient noise" looks to be below but the banging noise is above. In pseudocode, something like this:

raw_signal = importSignal("signal file path.wav")
signal_mean = mean(raw_signal)
signal_nSamples = length(raw_signal)
DC_removed = zeros(1,signal_nSamples)
for i = (0:signal_nSamples-1) {
    DC_removed[i] = raw_signal - signal_mean
}
signal_max = max(DC_removed)
threshold = 0.5*signal_max // or what ever you want it to be
clean_signal = zeros(1,signal_nSamples)
clean_signal_index = 0
// take out samples under the threshold
for i = (0:signal_nSamples-1) {
    if (DC_removed[i]<threshold){
        clean_signal[clean_signal_index] = DC_removed[i]
        clean_signal_index = clean_signal_index + 1
    }
}

Alternatively use a running mean (median or other filter, get creative if you like) and have the threshold be a level over this running mean. This is similar to above, but adapts to changes in background noise, so that the noise of a passing ship won't be filtered out by your threshold. If I had this file and it was short enough to do it all in the RAM I'd do a slightly more statistical approach:

  1. Do the running mean DC filter
  2. Find a suitable window length (e.g. a window length that gets a few bumps in each)
  3. For each window compute the mean of the squared samples (sample²)
  4. For each window find the value of e.g. 90th percentile, that is now your threshold
  5. Remove values above this threshold

Notes:

  1. You can get much more fancy that I was here when removing DC offset, e.g. use a running mean to subtract from the data to account for DC drift over the deployment.
  2. If your signal is too long to do this in the RAM of your machine, you have to get more creative and read the file piece-wise (python can do this, probably MatLab too). I think a tool like PAMGuard can detect impulses etc. and can handle very large files and exports detections to another file so you can use them to filter your signal.
  3. If you make a spectrogram of your sound you can use the spectral information in your bumps to trigger your threshold filter.

2022/09/15 Edit: I dug this octave/MatLab code out, I've used it in the past, works well enough to find start and end of impulses, you could use it to find the impulses and remove them.

function [start_samples, end_samples] = impulsefinder(signal, fs, step_in_sec, impulse_threshold_dB)
    pkg load signal

    sig_kurtosis = kurtosis(signal); 
    sig_is_impulsive = sig_kurtosis>3; # check if signal is impulsive, kurt>3 (gaussian noise kurt = 3)

    if sig_is_impulsive
        step_samples = ceil(step_in_sec*fs);
        impulse_threshold = 10^(impulse_threshold_dB/10); # make threshold linear
        
        # use running mean to identify start of impulses
        sig_hilbert = abs(hilbert(signal));
        sig_hilbert_running_mean = filter(ones(step_samples,1)/step_samples,1,sig_hilbert);
        clear signal # conserve memory
        clear sig_hilbert # conserve memory
        
        # make running mean and downsample
        step_to_down_sample = max(1,round(step_samples/10));
        mask = sig_hilbert_running_mean(1:step_to_down_sample:length(sig_hilbert_running_mean));
        clear sig_hilbert_running_mean # conserve memory
        mask_fs = fs/step_to_down_sample;
        
##        mask_time = (0:length(mask)-1)/mask_fs; # used i debug
##        plot(mask_time,mask) # used in debug
        mask_threshold = max(mask)/impulse_threshold
        exceedance_samples = find(mask >= mask_threshold);
        # go through exceedance_samples to find samples belongin to same impulse
        start_samples_mask(1) = exceedance_samples(1);
        for i=2:length(exceedance_samples)
            prev_sample = exceedance_samples(i-1);
            sample = exceedance_samples(i);
            n_starts = length(start_samples_mask);
            if ~(sample-1 == prev_sample)
                start_samples_mask(n_starts+1) = sample;
            else # this overwrites the curernt end_samples_mask a lot, but it's probably as quick as the alternative if checks
                end_samples_mask(n_starts) = sample;
            endif
        endfor
        ## convert mask samples to real-time samples
        start_samples = start_samples_mask.*step_to_down_sample;
        end_samples = end_samples_mask.*step_to_down_sample;
    ##    start_times = start_samples_real./fs
    ##    end_time = end_samples_real./fs
    else
        disp("signal is not impulsive, no impulses found.")
        start_samples = [];
        end_samples = [];
    endif
endfunction


[starts, ends] = impulsefinder(sig, fs, 0.05, 3.0); ## holds start and end samples for impulses
$\endgroup$
1
$\begingroup$

What I typically do in such situations is to simply mask the timeseries after detecting the unwanted signal with NaN. After that I would remove NaN from the data or use functions, where one can omitnan.

Concerning the detector, any detector would be appropriate (i.e, threshold detector)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.