How can I remove a pulse from a soundwave?
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Valentin Puiu
le 2 Mar 2023
Commenté : Image Analyst
le 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1312645/image.png)
Could you please tell me what am I missing here and how I could solve this?
Thank you.
1 commentaire
Réponse acceptée
William Rose
le 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 commentaires
William Rose
le 3 Mar 2023
Modifié(e) : William Rose
le 3 Mar 2023
[edit: Added figure showing the time domain signal without and with the added pulses.]
@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)');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1312825/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1312820/image.png)
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
le 3 Mar 2023
Modifié(e) : Walter Roberson
le 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?
Plus de réponses (2)
Star Strider
le 3 Mar 2023
Modifié(e) : Star Strider
le 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
le 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 commentaires
Image Analyst
le 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.
Voir également
Catégories
En savoir plus sur Get Started with Signal Processing Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!