# How can I remove a pulse from a soundwave?

1 view (last 30 days)

Show older comments

Valentin Puiu
on 2 Mar 2023

Commented: Image Analyst
on 13 Mar 2023

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

### Accepted Answer

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
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?

### More Answers (2)

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)

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

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).

.

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
on 13 Mar 2023

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!