Low Pass filter not working

I audioread() a signal and tried to apply low-pass filtering but it does not seem to have any change at all. The signal is a recording of lung sound and I wish to filter out the noise component.

[y,Fs] = audioread('mysound.wav')
Fs  = 44100;                                 % Sampling Frequency (Hz)
Fn  = Fs/2;                                 % Nyquist Frequency
Fco =   70;                                 % Passband (Cutoff) Frequency
Fsb =   100;                                 % Stopband Frequency
Rp  =    1;                                 % Passband Ripple (dB)
Rs  =   10;                                 % Stopband Ripple (dB)
[n,Wn]  = buttord(Fco/Fn, Fsb/Fn, Rp, Rs);  % Filter Order & Wco
[b,a]   = butter(n,Wn);                     % Lowpass Is Default Design
[sos,g] = tf2sos(b,a);
filt_sig = filtfilt(sos,g,y)

Neither the plot(), FFT() or soundsc() shows anything different. I've tried cheby filters as well. Am I doing anything wrong? Thanks for the help.

[Daniel M]

You would have to provide your data, or at least show the timeseries before and after, and the fft before and after filtering. As for soundsc, I suspect either your speakers do not play low-end frequencies very well, or the SNR of the signal is low below your cutoff frequency.

What does your noise consist of? This changes how you approach the filter.

Is it evenly distributed, white noise?

Is the noise separated from your signal spectrally?

Is the noise from a particular, characteristic external source?

[Yipp Chun Munn]

The frequency domain and Time domain. No matter what filter I apply, the Far right frequency components does'nt go away.

Answer 1

[Star Strider]

Your stopband attenuation is likely not sufficient to produce any difference.

Try this instead:

Rs  =   60;                                 % Stopband Ripple (dB)

Increase ‘Rs’ (to increase the stopband attenuation) as necessary to get the result you want.

Also, always use freqz (or similar functions) to see what your filter is doing:

figure
freqz(sos, 2^14, Fs)

[Yipp Chun Munn]

Hi, thanks for the input, I've tried playing with the Rs as well. There is no change as I gradually increase but then past a certain vaule suddenly the entire signal becomes 0.

[Star Strider]

My pleasure.

I have no idea why the entire signal should suddenly become 0 (the filter could have become unstable), however there are ways to produce filters that will only in extreme situations be unstable.

Try this:

Fs  = 44100;                                 % Sampling Frequency (Hz)
Fn  = Fs/2;                                 % Nyquist Frequency
Fco =   70;                                 % Passband (Cutoff) Frequency
Fsb =   100;                                 % Stopband Frequency
Rp  =    1;                                 % Passband Ripple (dB)
Rs  =   200;                                 % Stopband Ripple (dB)
[n,Wn]  = buttord(Fco/Fn, Fsb/Fn, Rp, Rs);  % Filter Order & Wco
[z,p,k] = butter(n,Wn);                     % Lowpass Is Default Design
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^14, Fs)
set(subplot(2,1,1), 'XLim',[0 500])         % ‘Zoom’ Frequency Range
set(subplot(2,1,2), 'XLim',[0 500])         % ‘Zoom’ Frequency Range

[Star Strider]

The entire code, including the Fourier transforms of the unfiltered signal, filter design, filtered signal, and the empirical transfer function:

[y,Fs] = audioread('mysound.wav');
Fs  = 44100;                                 % Sampling Frequency (Hz)
Fn  = Fs/2;                                 % Nyquist Frequency
Fco =   70;                                 % Passband (Cutoff) Frequency
Fsb =   100;                                 % Stopband Frequency
Rp  =    1;                                 % Passband Ripple (dB)
Rs  =   200;                                 % Stopband Ripple (dB)
[n,Wn]  = buttord(Fco/Fn, Fsb/Fn, Rp, Rs);  % Filter Order & Wco
[z,p,k] = butter(n,Wn);                     % Lowpass Is Default Design
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^14, Fs)
set(subplot(2,1,1), 'XLim',[0 500])         % ‘Zoom’ Frequency Range
set(subplot(2,1,2), 'XLim',[0 500])         % ‘Zoom’ Frequency Range

filt_sig = filtfilt(sos,g,y);

L = size(y,1);
Fn = Fs/2;
FTy = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTy(Iv))*2)
grid
title('Fourier Transform of Original Signal')

FT_filt_sig = fft(filt_sig)/L;
figure
plot(Fv, abs(FT_filt_sig(Iv))*2)
grid
title('Fourier Transform of Filtered Signal')

TF = FT_filt_sig ./ FTy;

figure
plot(Fv, abs(TF(Iv))*2)
grid
title('Transfer Function')

You may need to limit the frequency axes of the Fourier transform plot using xlim, for example:

xlim([0 250])

to see the region of the spectrum that demonstrates the filter result.

[Yipp Chun Munn]

First of all, very grateful for your help. After playing around with the passband and stopband frequency, I choose Fco = [50] and Fsb = [100] and got this result. I believe that is lung sound signal that I want(I researched and found it to be in the range of 60-600 hz). However when I soundsc(filt_sig, Fs) or sound(filt_sig, Fs) I'm not really quite able to hear anything but I believe something is playing. Could it have been my desired signal go attuenated? Thanks.

[Star Strider]

I always use sound, not soundsc, since soundsc scales the output to ±1. Using sound risks distorting the output, however it is possible to hear it. You can also multiply the argument matrix by a scalar (I usually use 10 to 100) to amplify it, then use sound (not soundsc) to hear it.

Otherwise, your filtered signal appears to have only one frequency at around 59 Hz. This is not a lung sound in my experience, since they cover a range of frequencies, and seems more like line-frequency noise. Since the 60 Hz hum could be overpowering everything else in your file, see Remove the 60 Hz Hum from a Signal to see if removing it allows you to hear the data in your file.

[Yipp Chun Munn]

Dear Star Strider, I've tried the Remove the 60 Hz hum from a signal solution but still to no avail. I've continued to try various filters still and the best possible outcome I've gotten is

This is was obtained using a cheby2 bandpass filter at Passband 850-900hz and stopband at 600 1000Hz. Not perfect at all but at least I can make out the 3 waveform above are the 3 breathing sounds. Also when I sound() this signal, the breathing is more audible. Here's the confusion, when I abs(fft()) the filtered signal, I get components at > 5khz. How should this be as it is clearly way out of my bandpass. Very sorry for asking so many question, thanks.

[Star Strider]

The last image you posted appears to show your filter working as it should. No filter will completely eliminate the frequencies outside the passband and stopband limits, however they will be significantly attenuated. That is what you want.

A better and more efficient filter solution would likely be an elliptic filter.

Prototype code for that would be:

Fs = 2250;
Fn = Fs/2;
Wp = [59 61]/Fn;                                                    % Stopband Frequency (Normalised)
Ws = [58 62]/Fn;                                                    % Passband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 90;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^18, Fs)                                                % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 100])                                 % Frequency Axis Limits
set(subplot(2,1,2), 'XLim',[0 100])                                 % Frequency Axis Limits

% signal_filt = filtfilt(sos, g, signal);                             % Filter Signal

Make appropriate changes for your sampling frequency, filter passbands and stopbands, and signal. Elliptic filters can handle very narrow transition regions, so:

Wp = [850 900]/Fn;                                                  % Stopband Frequency (Normalised)
Ws = [800 950]/Fn;                                                  % Passband Frequency (Normalised)

will produce a stable filter, likely with better charactersitics.

EDIT —

Corrected typographical error.

Answer 2

[Daniel M]

You show a figure that essentially looks like the following (also attached):

fs = 1000;
t = 0:1/fs:10;
x = sin(2*pi*10*t);
X = fft(x);
figure
plot(abs(X))

This is completely normal. It is the two sided amplitude spectrum.

You should read more about fftshift, and follow the example for the single sided amplitude spectrum on page for fft.

"No matter what filter I apply, the Far right frequency components does'nt go away."

That's because they're not supposed to.