How can I remove a pulse from a soundwave?

4 views (last 30 days)
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.
  1 Comment
Walter Roberson
Walter Roberson on 2 Mar 2023
Wouldn't you want a low-pass filter rather than a high-pass filter?

Sign in to comment.

Accepted Answer

William Rose
William Rose on 2 Mar 2023
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.
  3 Comments
Walter Roberson
Walter Roberson on 3 Mar 2023
Edited: Walter Roberson on 7 Mar 2023
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?

Sign in to comment.

More Answers (2)

Star Strider
Star Strider on 3 Mar 2023
Edited: Star Strider on 3 Mar 2023
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).
.
  1 Comment
Valentin Puiu
Valentin Puiu on 5 Mar 2023
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.

Sign in to comment.


Image Analyst
Image Analyst on 3 Mar 2023
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:
  3 Comments
Image Analyst
Image Analyst on 13 Mar 2023
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.

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!