How can I remove a pulse from a soundwave?

Hello,

I'm trying to remove a 10Hz pulse that I have added to a sound loaded from a file.

[x, Fs] = audioread('sample3s.wav');
t = (0:length(x)-1)/Fs;

pulse = max(square(2*pi*10*t, 10), 0);

noisy_x = x(:,1) + pulse.';

fstop = 10; 

order = 4;

[b, a] = butter(order, fstop/(Fs/2), 'high');

y_filtered = filter(b, a, noisy_x);
plot(t,y_filtered);

Despite using the filter, the pulses remain in the sound.

Could you please tell me what am I missing here and how I could solve this?

Thank you.

[Walter Roberson]

Wouldn't you want a low-pass filter rather than a high-pass filter?

Answer 1

[William Rose]

@Valentin Puiu,

Sound is in the range 20 Hz-20 kHz. A 10 Hz signal will be felt more than heard. You have added a square wave, which has power at 10, 20, 30, ... Hz, due to the usual harmonic series.

I agree with @Walter Roberson that the highpass filter is the way to go. That would preserve the frequencies in the audio range. But a highpass Butterworth filter with fco=10 Hz will only attenuate a sinusoid at 10 Hz by 3 dB, which is not very much attenuation. And the higher harmonics wil not be attenuated to any significant degree.

If you really want to remove the 10 Hz square wave, you could try a series of narrow notch filters at 10, 20, 30, 40... Hz. But narrow notch filters can do funny things. Caveat emptor.

[Walter Roberson]

Actually I suggested low-pass rather than high pass. When people talk about "pulses", I typically think of short-duration events like the numerous spikes in the signal diagram. I didn't pay close enough attention to the situation at hand.

[William Rose]

[edit: Added figure showing the time domain signal without and with the added pulses.]

Thank you , @Walter Roberson. I too did not pay close enough attention.

@Valentin Puiu, when I said "A 10 Hz signal will be felt more than heard", I should have said that a 10 Hz sinusoidal signal will be felt more than heard. The pulses you added are a square wave with 10% duty factor at 10 Hz, so they have lots of high frequencies to create the sharp transitions. The file will be the background sound plus brief clicks at a 10 Hz rate.

SInce you did not provide your sound file, I edited a sample to extract a C major triad, and added the pulse that you added. Refer to attached files.

[x,Fs]=audioread('cMaTriad.mp3');
t=(1:length(x))/Fs;
pulse = max(square(2*pi*10*t, 10), 0);
xN=x(:,1)+pulse';
yN=xN/max(xN);                      % to prevent clipping
audiowrite('cMaTriadN.mp4',yN,Fs);

Commands above read the file, add the pulses, and save the result.

subplot(211); plot(t,x,'-r'); grid on;
subplot(212); plot(t,yN,'-r'); grid on; xlabel('Time (s)');

The plots above show the signal before and after pulse is added.

Let's look at the spectrum of the pulses:

[pxx,f] = pwelch(pulse-mean(pulse),2*Fs,Fs,2*Fs,Fs);    % estimate power spectrum
figure; semilogy(f,pxx,'-b'); grid on; xlim([0 200])    % plot power spectrum
xlabel('Frequency (Hz)'); ylabel('Power'); title('Pulse Power Spectrum')

The code above makes the figure below:

The power spectrum shows that pulse has power at 10, 20, 30, 40, ... Hz, as I said it would, in my initial answer. You could try to get rid of this with a series of notch filters, but I think the results will not be satisfactory. I think that because notch filters introduce ringing and overshoot, and because you would probably need at least 5 notch filters to significantly attenuate the 10 Hz clicks, which will exacerbate the ringing and overshoot problems.

[Walter Roberson]

comb filter might make the code easier, https://www.mathworks.com/help/dsp/ref/fdesign.comb.html but will not necessarily help.

If you were working entirely digitally I wonder if there would be a useable way to figure out the amplitude and duty cycle of the square wave and then synthesize a wave and subtract it out?

Answer 2

[Star Strider]

If the pulse is at a single frequency with known bandwidth, just use a bandstop filter. If there is energy at several frequencies, a comb filter would work. These are relatively easy to code using kaiserord and fir1, the disadvantage being that FIR filters are long filters and require long signals. Use the filtfilt function to do the actual filtering.

EDIT — (3 Mar 2023 at 15:13)

Added simulated original signal, pulses, and filters —

This is likely going to be a difficult signal to work with, regardless —

Fs = 250;
Ls = 3.25;
t = linspace(0, Ls*Fs, Ls*Fs+1).'/Fs;
x = sin(2*pi*5*t)/5;
pulse = max(square(2*pi*10*t, 10), 0);
noisy_x = x(:,1) + pulse;

figure
plot(t,noisy_x)
xlabel('t')
ylabel('Amplitude')
title('Original Signal')

Fn = Fs/2;
L = numel(t);
NFFT = 2^nextpow2(L)

NFFT = 1024

FTnx = fft((noisy_x-mean(noisy_x)).*hann(L), NFFT);
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTnx(Iv))*2)
grid

[pks,locs] = findpeaks(abs(FTnx(Iv))*2, 'MinPeakProminence',1);
hold on
plot(Fv(locs), pks, '^r')
hold off
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform of Original Signal')

% return

fcomb = reshape([[-5 -3 3 5]*0.05+Fv(locs).'].', 1, []);                % Design FIR Comb Filter
mags = [[1 0 1] repmat([0 1],1,numel(locs)-2), [0 1]];
dev = [[0.5 0.1 0.5] repmat([0.1 0.5],1,numel(locs)-2), [0.1 0.5]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
n

n = 2100

hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
 
figure
freqz(hh,1,2^20,Fs)

% return

nxf = repmat(noisy_x, ceil(n*3/numel(noisy_x)), 1);                     % Lengthen Signal
nx_filt = filtfilt(hh, 1, nxf);                                         % Apply FIR Comb Filter
nx_filt = nx_filt(1:numel(noisy_x));                                    % Shorten Signal To Original Length
nx_filt1 = lowpass(nx_filt, Fv(locs(1))*1.1, Fs);                       % Apply 'lowpass' Filter To Comb-Filtered Signal
nx_filt2 = lowpass(noisy_x, Fv(locs(1))*1.1, Fs);                       % Apply 'lowpass' Filter To Original Signal

figure
plot(t,noisy_x)
hold on
plot(t, nx_filt1)
hold off
grid
xlabel('t')
ylabel('Amplitude')
title('Comb-Filtered & Lowpass-Filtered Signal')

figure
plot(t,noisy_x)
hold on
plot(t, nx_filt2)
hold off
grid
xlabel('t')
ylabel('Amplitude')
title('Lowpass-Filtered  Signal')

It may not be possible to completely recover the original signal (without the additional pulses).

.

[Valentin Puiu]

Hello @Star Strider,

Unfortunatelly I can accept a single answer to the question, but I would like to thank you for your effort and for the great insight that you gave me.

Thank you.

Answer 3

[Image Analyst]

Well using a frequency filter as the others suggested is not what I'd try with pulses in your signal that looks like that. It looks like any part of the signal that is above 0.5 or below -0.5 is what you want to remove, so I'd simply filter it in the time domain with thresholding and masking:

mask = abs(noisy_x) > 0.5;
repaired_x = noisy_x; % Initialize.
% Now remove elements
repaired_x(mask) = []; % or = 0 if you want to leave in those samples but just silence them.

If you have any more questions, then attach your data, 'sample3s.wav' with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Image Analyst]

Not sure why this didn't work, but I guess we'll never know because you won't attach your data.

[Valentin Puiu]

Hello,

Unfortunately, I can mark a single answer as solution, even if more qualify.

When I made the question I didn't even tought that an approach like the one you proposed is possible, otherwise I would have specified that I have to remove the pulse using a FIR filter with window method. Also, I was not looking for a 100% solution (the entire problem - school assignment - is slightly more complex), but to a guidance in how to solve it. The answers from @William Rose and @Star Strider taught me that the pulse has multiple harmonics and that I should filter all of them.

The data didn't mattered. The file was just a random sound file, short enough to allow for a fast simulation and the pulse was added on it. Thus, I considered that any sound file is enough.

Thank you for taking your time to provide an answer and I hope that this answer sheds some light to why another answer was prefered.

[Image Analyst]

OK, just as long as you know that filtering or smoothing out the spikes (like the others did) is not the same as removing them entirely (like I did). The signals will look different.