Eliminating noise using butterworth bandpass filter

The aim of this task is to build a bandpass filter to filter out the noise from a given data. The sampling frequency was 3490Hz. I have used fft to determine what the noises are in the given signal which are: 20Hz. 110Hz, 200Hz and 470Hz. The peak at 20Hz is peculiar to me though as the left hand side of the fft figure should be a mirror image of the right hand side of the figure. This is the fft result that I achieved:

The frequency is in bins and was converted to Hz by comparing with the sampling frequency.

Next, I need to determine the actual noise of the signals. I am not sure how to determine this. After all that, I can then possibly make a bandpass filter by the basic code flow below:

Fs = 3490;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [futcutlow fcuthigh]/Fn;                               % Passband Frequency (Normalised)
Ws = [futcutlow*0.95 fcuthigh/0.95]/Fn;                     % Stopband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs = 150;                                                   % Stopband Ripple (dB)
[n,Wn] = buttord(Wp,Ws,Rp,Rs);                              % Filter Order
[z,p,k] = butter(n,Wn);                                     % Filter Design
[sosbp,gbp] = zp2sos(z,p,k);                                % Convert To Second-Order-Section For Stability
freqz(sosbp, 2^20, Fs)                                      % Filter Bode Plot
filtered_signal = filtfilt(sosbp, gbp, signal);             % Filter Signal

The problems that I've encountered is: 1) Why is there a spike in the left hand side of the fft at 20Hz? 2) How do I determine the actual noise that I need to eliminate? 3) How can I make the basic code flow work from that point onwards?

Answer 1

[Star Strider]

Assuming ‘signal’ is a vector (as it appears to be in the image), you can accurately calculate the fft (link) with:

L = length(signal);                                          % Signal Length
FTsignal = fft(signal)/L;                                    % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                          % Frequency Vector
Iv = 1:numel(Fv);                                            % Index Vector

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

That information will produce the correct amplitude-frequency data, and will allow you to design your filter.

For your signal, a Chebyshev Type II filter may work better, if you want to eliminate or pass specific frequencies. The filter design then becomes:

Wp = [fcutlow fcuthigh]/Fn;                                 % Passband Frequency (Normalised)
Ws = [fcutlow-1 fcuthigh+1]/Fn;                             % Stopband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws]  = cheb2ord(Wp,Ws,Rp,Rs);                            % Filter Order
[z,p,k] = cheby2(n,Rs,Ws);                                  % Filter Design, Sepcify Bandpass
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability

figure(3)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot

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

You may need to experiment to get the result you want.

[Hafiz Ahmad]

signal = importdata('A-Scan.txt');
Fs = 3490;
Fn = Fs/2;
L = length(signal);                                          % Signal Length
FTsignal = fft(signal)/L;                                    % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                          % Frequency Vector
Iv = 1:numel(Fv);                                            % Index Vector
figure(2)
plot(Fv, abs(FTsignal(Iv))*2)
grid

This is the graph that I attained by using the code above. However, I still don't quite understand how this fft is used to design the filter.

[Star Strider]

You design the filter based on the results the fft provides. For example, there is a peak at about 100 Hz, so if you want to specifically pass it and reject the other data, the filter passband could be 75 Hz to 125 Hz, depending on what information you want.

Having the correct frequency-amplitude information is important for correct filter design.

[Hafiz Ahmad]

How do I input the signal into the filter?

[Star Strider]

Use the filtfilt function:

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

[Hafiz Ahmad]

signal = importdata('charlie.txt');
fcutlow = 160;
fcuthigh = 320;
Fs = 3490;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [fcutlow fcuthigh]/Fn;                               % Passband Frequency (Normalised)
Ws = [fcutlow*0.95 fcuthigh/0.95]/Fn;                     % Stopband Frequency (Normalised)
Rp = 160;                                                   % Passband Ripple (dB)
Rs = 320;                                                   % Stopband Ripple (dB)
[n,Wn] = buttord(Wp,Ws,Rp,Rs);                              % Filter Order
[z,p,k] = butter(n,Wn);                                     % Filter Design
[sosbp,gbp] = zp2sos(z,p,k);                                % Convert To Second-Order-Section For Stability
freqz(sosbp, 2^20, Fs)                                      % Filter Bode Plot
filtered_signal = filtfilt(sosbp, gbp, signal);             % Filter Signal
imwrite(uint16(filtered_signal - double(intmin('int16'))), 'x4.png');

I am trying to run it using butterworth filter this time but I can't simply run it without going into multiple errors. Can you help me where the error lies?

I am storing it as .png because it is an image.

[Star Strider]

The filter design runs without error, although it produces a very long filter, with ‘n=130’. Consider instead using:

Ws = [fcutlow*0.90 fcuthigh/0.90]/Fn;                     % Stopband Frequency (Normalised)

If you need a shorter transition region, use a Chebyshev filter.

I have no idea what ‘charlie.txt’ contains. If you are filtering an image, you must define the filter differently, use the fspecial function to create it as a 2D filter, and use the imfilter function to filter it.

I do not usually do much image processing, so I do not use these functions much, and have not used them recently. Note that the fwind1, ftrans2, and other image filters require FIR filters, not the IIR filter you are designing here.

You will benefit by searching the File Exchange for demos Image Analyst has posted. (I cannot find any just now that specifically address filtering.) Also see Steve on Image Processing (link) for information on how to do image filtering, and the examples in Image Filtering and Enhancement (link).

[Hafiz Ahmad]

charlie.txt This is the file that I want to filter. Using Chebyshev filter, I can get the desired filtered image (I'm not sure how). Now, I'm trying to compare using butterworth.. which gives an error of not enough samples.

[Star Strider]

When I downloaded your (349x4007) double file and then tried to visualize it with image and imshow, it produced nothing but garbage. With mesh, it seems to have an initial oscillation and some peaks, and then decays to an equilibrium condition at a level of about 5000. I did not attempt to filter it.

I suspect the reason the Butterworth design failed is that there are not enough rows for such a long filter, since filtfilt operates along the first array dimension of ‘signal’ with size greater than 1. (A Chebyshev design with those passbands and stopbands will produce a shorter filter.) Transpose ‘signal’ to be (4007x349) and try again, with both filters.

[Hafiz Ahmad]

Hello again! could you explain what's the difference between frequency and ripple? What value is input into the ripple (passband ripple and stop band ripple)?

[Star Strider]

The frequencies define the passbands and stopbands.

Chebyshev filters have either passband ripple (Type I) or stopband ripple (Type II). The ripple is a small variation of the passband or stopband. In the stopband, the ripple also defines the attenuation. (In Butterworth designs, the stopband ‘ripple’ defines the stopband attenuation, since Butterworth filters have maximally flat passbands and stopbands.) Consult any good textbook on signal processing for more extensive discussions of filter design.

[Hafiz Ahmad]

[n,Ws]  = cheb2ord(Wp,Ws,Rp,Rs);

Here, is Rs the stopband attenuation? I have read from design filter that it is stopband attenuation. From an A-scan view of the data as can be seen below, the attenuation would be from the highest peak to the highest point of the ripple at the negative part of the signal, by calculation that would be 2.8x10^4 + 2.15x10^4 = 49500dB. This value gives an error. --- (highest peak to the highest point of the ripple) So my questions are: (1) Is Rs the attenuation? (2) A flat Rp is assumed to be 1? (3) How exactly do we find Rs?

[Star Strider]

In your Comment, you are analysing the filtered signal, not the filter characteristics. You have to compare the filtered signal to the unfiltered signal (preferably using the Fourier transform of each) to understand what the filter is doing.

To see the filter characteristics, use a Bode plot, given by the freqz function for discrete filters. That will tell you what the filter is doing from d-c to the Nyquist frequency, since it compares the Fourier transforms of the input and output of the filter at each frequency.

You do not ‘find’ ‘Rs’ or any other filter design parameter (although you can calculate the optimal order based on other design parameters). You specify ‘Rs’ in the design of your filter, depending on what you want the filter to do. If you want the filter to have a 50 dB stopband attenuation, specify:

Rs = 50;

in your design code.

[Hafiz Ahmad]

How to we determine the Rp and Rs then? Do we just trial and error? or is there other ways to determine it?

[Star Strider]

They are design parameters, and so are determined by what you want the filter to do, and your constraints.

They are more important in Chebyshev and other such designs. In Butterworth designs, they simply determine the passband and stopband attenuation. Their values are determined in hardware designs by the complexity of the filter circuitry and the available components, and in software designs by filter length and computational speed considerations. If none of those are important, the values can be anything you want, although in practical terms ‘Rs’ usually only needs to be about 50 dB.

[Hafiz Ahmad]

When I use fcutlow to be 322 and fcuthigh to be 563, it gives this:

Error using filtfilt>getCoeffsAndInitialConditions (line 157) Data length must be larger than 366 samples.

Error in filtfilt (line 97) [b,a,zi,nfact,L] = getCoeffsAndInitialConditions(b,a,Npts);

Error in Bscanfilter (line 13) signal_Filtered = filtfilt(sos, g, signal);

What is the reason for this?

[Star Strider]

You have a long filter, and apparently a short signal.

The frequencies may be very close together with respect to ‘Fn’, in the context of your other design parameters.

Experiment with other tweaks to your filter to get it to work with your signal, and do what you want it to do.

[Hafiz Ahmad]

Why is Chebyshev Type II more preferable than butterworth and Chebyshev type I? I've read about the frequency dependance of Chebyshev filters but I am still unclear on why type II is preferred. Thanks!

[Star Strider]

Like all filter choices, it depends on the intended application. A Chebyshev Type II filter has a flat passband, a narrow transition region, and ripples in the stopband. (A Chebyshev Type I has a flat stopband and ripples in the passband.) A Butterworth filter has maximally flat passband and stopband, and a more gradual transition region.

I usually prefer a Chebyshev Type II filter because signals of biomedical processes are characteristically low-frequency band-limited (typically 0-500 Hz), so usually require narrow transition regions for adequate signal separation. For other signal processing applications (for example audio, where frequencies are in the range of 0-20000 Hz, and where signal distortion can be more noticeable), a Butterworth filter is likely a better choice.

Filters implemented in hardware have different options because they have different constraints. For example, a Bessel design is best for a hardware anti-aliasing filter and some other applications, since it has a linear phase characteristic, and so does not distort the signals in its passband. Other designs cause phase distortion.

[Hafiz Ahmad]

What about when compared to chebyshev type I? What is the decisive factor to pick type II?

[Star Strider]

If you want a flat stopband and can tolerate ripples in the passband (usually designed to be about 1 dB), choose Type I. If you want a flat passband and can tolerate ripples in the stopband, choose Type II.

[Hafiz Ahmad]

I have 3 more questions to ask: (1) Why is filter order 2 more stable than higher orders? (2) When is elliptic filter favorable to be used in a filter design? (3) Higher frequencies are more readily absorbed or scattered than lower frequencies.. is this correct? Thanks for your help so far.

[Star Strider]

(1) The filter design and implementation governs the stability, not the order, although more complicated filters can be more difficult to implement. It is probably difficult to design and implement an unstable order 2 filter.

(2) Elliptic filters produce the smallest order filter for specific design parameters, and have the smallest transition bandwidth. If such efficiency is important, an elliptic filter would be the better option.

(3) I don’t understand the question. However I doubt that there is a difference.

[Hafiz Ahmad]

Why is there an initial sudden peak in fourier transform? but when I don't use the Nyquist frequency but the sampling frequency, the peak is not present. Why is this?

[Star Strider]

I have no idea what you’re referring to. My guess is that you’re seeing the d-c (constant) offset, that appears at 0 Hz. You can eliminate it by taking the fft of your signal after subtracting the mean of your signal:

FTsignal = fft(signal - mean(signal))/length(signal);

Always use the Nyquist frequency. It is the highest frequency at which you can uniquely resolve a signal without aliasing.