Implementation of Narrow band pass filter ( Butterworth)

Hi all,

I'm struggling with implementation of narrow band pass filter. I understand from few suggestions (feedback) that it's too narrow and impulse is becoming too large. I also tried using the decimation but even it didn't work. I've posted my question in stackoverflow. Here is the link to it.

http://stackoverflow.com/questions/29111768/band-pass-implementation-matlab?noredirect=1#comment46453099_29111768

Could someone help me to solve this. Since this is part of my project. I'm little running short of time.

Thanks in advance for all your help and suggestion.

Answer 1

[Star Strider]

This seems to work, although freqz has problems with it for some reason:

% DEFINE FILTER PARAMETERS 
Fs  = 2E+6;
Fn  = Fs/2;
Fc1 =  630;
Fc2 =  640;
Fs1 = Fc1*0.8;
Fs2 = Fc2/0.8;
Rp  =    1;
Rs  =   10;

% DESIGN FILTER
[n,Wn]  = buttord([Fc1 Fc2]/Fn, [Fs1 Fs2]/Fn, Rp, Rs);
[b,a]   = butter(n,Wn);
[sos,g] = tf2sos(b,a);

% CHECK FILTER PERFORMANCE
% figure(1)
% freqz(sos, 512, Fn)
s = zeros(1, Fs);
s(fix(Fs/2)) = 1;                       % Create Impulse Signal
y = filtfilt(sos, g, s);                % Filter Signal
tfs = fft(y)./fft(s);                   % Calculate Transfer Funciton
Fv = linspace(0,1,length(s)/2+1)*Fn;    % Frequency Vector
ix = 1:length(Fv);                      % Index Vector

figure(2)
plot(Fv, abs(tfs(ix)))
grid
axis([0  1000    ylim])

[Bharath]

Thanks for your answer. I tried to run this code but still I've unclear orbit plots. Also could you please explain me why did you use Fc1*8 and Fc2/0.8 and how did you calculate Fs1 and Fs2?

In the end the filter coeffcients are converted by using tf2sos into some variable sos and g. Can they be interpreted in filtering the signal with 'filter' command too or only 'filtfilt' has to be used?

Also I didn't undestand your last set of code. Was it to just show the filtered and unfiltered signal difference?

tfs = (fft(y)./fft(s))/length(s);
Fv = linspace(0,1,length(s)/2+1)*Fn;
ix = 1:length(Fv);

Also I'm attaching the orbit plots using this set of code.

As you can see there is lot of other frequency components which are not giving clean orbit.

Thanks for your time and answer.

[Star Strider]

My pleasure.

I leave the orbit plots to you. I’m just sorting out your filter.

There’s a certain ‘art’ to filter design. As a first approximation, I calculate the stopband frequencies (‘Fs1’, ‘Fs2’) for any filter as 0.8 times the low passband and 1.25 times the upper passband. (This usually works; if it doesn’t, I adjust them.) The buttord function uses these to determine the appropriate filter order.

I found significant problems with your filter design, and I gave you code to design and implement it efficiently to get the result you want. I converted your filter to a second-order-section representation (with tf2sos) for stability. This is common practice, and the butter function discusses it in detail.

Since you have the Signal Processing Toolbox, use filtfilt. It does not introduce any phase distortion, while filter does. (All this is explained clearly in the documentation for the various functions, so I won’t go into it in detail here.)

The three ‘tfs’, ‘Fv’, and ‘ix’ assignments calculate the frequency-domain transfer function, the frequency vector for the plot, and the index vector for the plot. (The index vector just makes the plot call easier to write.) Because your passband is very small in relation to your sampling frequency, the freqz function did not show the passband with sufficient resolution, so I created my own version. (Note: I edited my original code to provide a much better ‘s’ and transfer function calcualtion.)

[Bharath]

Thanks a ton! "That's true filter design is an art and you are a masterclass artist." Your answer was very informative. I just need to clarify few other points. kindly bear with me.

1. Once I've [sos,g] values calculated....all I need to do is filter my signal using the

y = filtfilt(sos, g, s); % where s is my raw signal?

• Is the impulse response and transfer function only used to plot the filter response right?
• How are the Rp and Rs values varied?
• I also tried varying the values for 'Fs1' and 'Fs2' sometimes it gives me an error saying they can't have overlapped values.
• After applying this filter to my raw signals X and Y. I seem to get distorted shapes. Although the filter impulse response shows its working fine. What do you think the problem could be with? I'm attaching the image for your reference. As you can see I still don't understand if the filter is really doing the job.

• Few people commented saying my sampling is a over kill for the frequency I'm filtering into. what is your suggestion about this?

Sorry for many questions. I need to clear my mind when I wanted to learn something in detail.

Thanks for your answers and patience.

[Star Strider]

My pleasure! I very much appreciate the compliment.

1. Correct. My ‘s’ variable is your raw signal.

2. Is the impulse response and transfer function only used to plot the filter response right?

Correct. I always test my filter designs first to be sure the result is close to what I want.

3. How are the Rp and Rs values varied?

They can be any reasonable values. They are really not applicable to a Butterworth design, because it has relatively flat passband and stopband characteristics, but buttord wants them. I always choose 1 dB for the passband ripple, but the stopband ripple can be much higher. I choose 10 dB, but 30 dB is reasonable.

4. I also tried varying the values for 'Fs1' and 'Fs2' sometimes it gives me an error saying they can't have overlapped values.

The ‘Fs1’ value always has to be strictly less than the ‘Fs2’ value. If they are too close together, the filter becomes unstable, and even a second-order-section representation will not produce a usable filter. I would narrow the stopband frequencies instead, if you want a narrower passband. Using a factor of 0.9 produces a stable filter with a much narrower passband:

Fs1 = Fc1*0.9;
Fs2 = Fc2/0.9;

See if that improves your filtered signal.

5. The problem is likely that the filter passband is too wide. See #4 for my suggestion on narrowing it. Experiment with higher numbers as well, perhaps up to 0.95, until you get a filter that is stable and that also gives you a filtered signal that is much closer to what you want. If you have a signal that is very close to the one you want and you cannot successfully filter it with the passband filter, consider notching it out with cascaded narrowband bandstop (‘notch’) filter. You would design it the same way as your bandpass filter, except to specify 'stop' as your filter type. (See the documentation for Bandstop Butterworth Filter for details.) The two-filter implementation for this would be to filter your raw signal with your bandpass filter first, then filter the output of that filter procedure with your bandstop filter to get your final filtered signal. You can also experiment with a Chebyshev Type II filter design. It has much steeper stopband attenuation than a Butterworth filter.

Don’t apologise for asking questions! That is what MATLAB Answers is for. I just have to have the time to answer them in as much detail as I feel they need. Right now, I do.

[Bharath]

Thanks again for your answers.

1. I tried to experiment with different variables mentioned in your #4 point but that doesn't seem to filter more further.
2. I was trying to explore your #5 point 'Usage of bandstop filter' and ''Cheby2' filter. I tried using the bandstop filter twice to filterout frequencies below and above the frequency of interest.
3. They seem to filter and I'm getting little more clear orbits. But how do I know my filters are doing the right job. Should I test it using 'freqz' command and see how the plot is?
4. And also if I apply many filters does that decrease my amplitudes of signal?

5. I was using hanning window at the end of filters and that too looks like doing little job. In overall I've been using many filters. Is that a good approach in signal processing?

[Star Strider]

My pleasure.

Your first two items don’t seem to require a reply.

3. Unfortunately, freqz has problems with your filters because of the high sampling frequency and low passband frequency. Otherwise I would suggest it. Using my version is probably better for your application. (In that instance, ‘y’ is the last output of all your filter operations. The rest of my transfer function calculation and plotting code is unchanged.)

4. It could, but that is simply a problem in any signal processing application. If you want a ‘clean’ signal with an extremely narrow passband, you have to live with such tradeoffs.

5. Generally, using the fewest filters possible is the best solution. Usually you can do everything you need to with at most two or three cascaded filters. (I generally don’t use windows with filters because they usually don’t add anything.) One possible approach to produce an extremely narrowband filter is to cascade two Chebyshev Type II filters — one lowpass and one highpass — and adjust the overlap. This is not the ideal solution, but it may be the only workable option if you have an extremely narrow band of frequencies you want to filter. It should be stable, where an extremely narrow bandpass filter might not be.

[Bharath]

In the end I would say " Your answers and suggestions were more understandable than those bulky theoretical textbooks" Finally, I got a feeling that filter design variables and orders has to be something like 'trial & error' The moment we feel that's right then that's right...that comes more from experience and experimenting. Thanks a lot for the enlightenment . Meet you soon again with another question ;)

[Star Strider]

Thank you very much! The textbooks and reference books are intended to cover all potential applications, giving theory as well as practical examples. I still use them often for MATLAB Answers questions. You had already designed your filter, but needed help in getting it to work correctly, so that made it easier for me.

A lot of engineering is ‘trial and error’, the reason it is always necessary to ‘breadboard’ designs first, then tweak them until they do what you want.

I’ll do my best to answer it!

[Bharath]

I forgot to ask about this error which I get when using filtfilt command.

Warning: The first two inputs are interpreted as an SOS matrix and G vector. Use column vectors if you meant to specify a
transfer function. 
> In filtfilt>getCoeffsAndInitialConditions (line 126)
  In filtfilt (line 97)
  In Test_run_chin_half_spectrum_2 (line 61) 
Warning: The first two inputs are interpreted as an SOS matrix and G vector. Use column vectors if you meant to specify
a transfer function. 
> In filtfilt>getCoeffsAndInitialConditions (line 126)
  In filtfilt (line 97)
  In Test_run_chin_half_spectrum_2 (line 62) 
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  5.552883e-17. 
> In filtfilt>getCoeffsAndInitialConditions (line 200)
  In filtfilt (line 97)
  In Test_run_chin_half_spectrum_2 (line 67) 
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  5.552883e-17. 
> In filtfilt>getCoeffsAndInitialConditions (line 200)
  In filtfilt (line 97)
  In Test_run_chin_half_spectrum_2 (line 69) 
>>

Is this something to be ignored? Also in other context like butterworth pass filter can we use filtfilt instead of filter command?

[Star Strider]

Ignore it.

I got the same warning each time, and a perfectly understandable plot. I checked the filtfilt documentation and my code, and couldn’t find any problems. It may be due to the ‘SOS’ matrix reflecting the narrow bandwidth (10 Hz), low centre frequency (635 Hz), and high sampling rate (2E+6 Hz).

It’s a warning, not an error, so accept it as such. You’re asking the filter design functions and filtfilt to do some calculations that are likely close to the double-precision limits of MATLAB. That’s essentially what it’s telling you.

[Bharath]

Alright,Thanks!

[Star Strider]

My pleasure!

This is an interesting problem for me, too!

[Bharath]

Hi again! When I'm experimenting with 'cheby2' filter. It gives me good results but everytime I need to vary these parameter 'n'(Order) and 'Rs' (Stopband Attentuation). So that is becoming like many iterative process. Isn't there any kind of thumb rule for defining these parameters? I'm using this filter for stop band to cut out the unwanted frequencies above and below the interested frequency.

 [z,p,k] = cheby2(n,Rs,[1/(fs/2),630/(fs/2)],'stop');
 [sos,g] = zp2sos(z,p,k);
 X = filtfilt(sos,g,X);
 Y = filtfilt(sos,g,Y);

[Star Strider]

Hi!

I thought I mentioned this, but apparently I didn’t. Use the cheb2ord function to determine the order and passband/stopband for each of the Chebyshev filters. The Chebyshev filter design differs a little from the Butterworth design. You’re correct in using the SOS implementation.

The filtfilt implementation remains unchanged.

If I remember correctly, there used to be a section on each type of filter, discussing the design of each. I can’t find it in the new documentation.

[Bharath]

Is it something like this? I want to cut off frequencies between 10 - 780 Hz. But again Rp and Rs values we need to vary according to our needs? (trail and error?)

 % Define Bandstop Filter Parameters.
 Wp = 10/(fs/2);
 Ws = 780/(fs/2);
 Rp = 1;
 Rs = 10;
 [n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);
 [z1,p1,k1] = cheby2(n,Rs,[1/(fs/2),750/(fs/2)],'stop');
 [sos2,g2] = zp2sos(z1,p1,k1);
dataInX = X; 
X = filtfilt(sos2,g2,dataInX); %filter command filters 
dataInY = Y;
Y = filtfilt(sos2,g2,dataInY);

[Star Strider]

That doesn’t look like a bandpass or bandstop design to me. Just like the Butterworth, a bandpass/bandstop design requires a two-element vector for both ‘Wp’ and ‘Ws’, calculated the same way as for the Butterworth design. See the cheb2ord documentation for details.

Trial and error in real-world situations predominates, but you need to be designing the filter you intend! If you want a bandpass or bandstop filter, you need two-element ‘Wp’ and ‘Ws’. If you want to cascade and overlap a highpass and lowpass to get a narrower passband, design them separately (only one ‘Wp’ and ‘Ws’ value for each).

The cheb2ord and cheby functions are essentially analogous to their respective Butterworth functions.

[Bharath]

Thanks! It works now. Again going back to our previous discussion about amplitudes. I've got now a cleaner Orbit.But as you can see in the figure the axis scale has been tremendously increased beyond imagination.

Unfiltered Orbit

Filtered Orbit ( Please note the axis scale which is 1E+58.

Why does this behavior happen. As I increase my speed it kept on increasing something like 1E+250 too. ( Isn't that too large)? Is it something to do with the filters?

[Star Strider]

I can’t comment on the amplitudes, since I’m not certain how you’re acquiring or processing your signal, other than filtering it. (One observation: you might want to use axis equal to be certain the same axis increments are used on both.) It is quite possible that as you increase the speed, the amplitude increases. I would look at your raw signal for that information. (The easiest way to do that is to calculate the root-mean-square amplitude of the entire signal.)

I doubt it’s the filters, unless the transfer function plot maximum is greater than 1 (it wasn’t in my filter calculations). I divide the fft of the output by the fft of the input, so any amplification or attenuation should show up in that plot. It doesn’t, or at least it didn’t in my calculations.

[Bharath]

Thanks for your answer! I'm getting raw data set from external client. Probably, I guess they are recording it with a Yokagawa DL560 Scopecorder. I#m processing the signal in the following way.

1. Reading the raw signal from Data sets with sampling frequency defined in the Data sets.

2. Detrending the Raw signal to remove the DC offset.

3. Apply Low pass filter ( Butterworth filter) to remove the unwanted higher unwanted frequency components.

4. Apply Bandpass filter (Butterworth filter) to pick up the interested frequency component with certain bandwidth.

5. Apply Bandstop filter ( cheby2 filter) to remove the other frequency components above and below the interested freq.

6. Apply hanning window.

7. Plot the orbits for X and Y raw signals.

[Star Strider]

I’m still guessing that the amplitude of your raw signal is the problem.

If your raw signal is x (after detrending to remove the D-C offset), take the RMS of it with:

RMS = @(x) sqrt(mean(x.^2));

That should give you some idea of its magnitude. (For example, the RMS value of a sine wave with peak amplitude 1 is 0.707.)

[Bharath]

I checked the RMS value of my signal and it is

rms(X_detrend)

ans =

     3.818678867784307e-02

>> rms(Y_detrend)

ans =

     5.338664908137016e-02

What also I noticed is when I apply one filter say one band pass filter its fine... But when I try to apply also cascade of low pass and high pass filter using cheby2 filter it goes very high.

Apart from that everything seems fine and when the scales are weird. I'm little confused in interpreting the results. Thanks in advance for your help

[Star Strider]

To the best of my knowledge, passive filters shouldn’t amplify. (These discrete filter realisations are derived from passive filter prototypes.)

What do your filter cascades do with my ‘CHECK FILTER PERFORMANCE’ code? (You will have to modify it slightly to include your filter cascades so that ‘s’ is the input, and ‘y’ the output of the cascaded filters.) That code divides the output by the input, so any amplification should show up there.

[Star Strider]

Bharath’s Answer moved here ...

Thanks for your reply. I tried testing using the "filter performance check" code. I don't see any change in the scales of the axis. Also, I checked my raw signals and started plotting the results after applying the filters sequentially. What I figured was after the application of the filters the scales seems to change that too with the higher filter orders and If I don't use the higher filter orders I don't get the filtered orbits. I' understanding where I'm going wrong. Is the way I'm applying the filters wrong or something wrong with my fundamental signal processing.

[Star Strider]

My pleasure.

I’m not certain what the problem is. You might apply them in a different order, using two overlapping Chebyshev filters (high- and low-pass) instead of the bandpass and bandstop filters, or perhaps designing hardware filters to filter the signal before you digitise it. It depends on what you want to do.

If the ‘FILTER PERFORMANCE CHECK’ code did not reveal any amplification, I’m at a loss to explain it. If it gives you the spectrum response you want, a 1.4x amplification may be acceptable if it applies to all your signals. Consult with signal processing experts as your university. They may have ideas or solutions I haven’t thought of. Another possibility is to take the RMS values of your entire input signal and normalise your output signal to that.

[Bharath]

Hi again. Why is there a change in the signal pattern after applying filters. Aren't the filtered signals supposed to be pure sinusodial waves after filtering?

As you can see in the figure there has been change in amplitude for filtered signals. Why does this behavior occurs?

[Star Strider]

I don’t know what your filter inputs are or how you designed your filters, so I can’t say what the outputs should be. Very narrow passbands or very low frequency passbands can produce such results.

I would do an fft on the inputs and outputs to verify that your filters are doing what you want them to do.

[Bharath]

Thanks for your answer. My filter is designed according to code you shared with me some days back. Also, I tried checking the filter response using the freq function program you shared with me. I'm sorry... I should have mentioned this in my question.

[Star Strider]

The filter I designed works, or I’d not have posted it. My only conclusion is that it is not appropriate for the signal you are processing with it. The sampling frequency has to match, and there has to be a frequency component in your signal for the filter to work with.

[Bharath]

Thanks again for your reply. The signals I'm using are the same signal...I used previously for this filter application. Also, I've shared my data in my question too. Its the same data set... plotted in time domain with filters applied to extract one frequency and same data set to extract another freqeuncy. I used the bandpass butter filter.

[Star Strider]

‘...same data set to extract another freqeuncy.’

That’s likely the problem. You have to design and test an entirely new filter to extract the second frequency. Filter design is heuristic and can be challenging, especially with the narrow passband filters you seem to need. You can use the essence of my code to help you with your design, but you have to determine the passband and stability of the second-order-section (SOS) implementation before you use it.

It is important to note that your wanting to design a filter to do a particular task does not mean that design is feasible. It might actually be easier to design and implement your filters in hardware and recording those filtered signals separately (along with your unfiltered signal) than to do the processing digitally. Narrowband hardware filters are easier to design and implement (at least in my experience) than narrowband digital filters.

[Star Strider]

Harry — Were I to design that same filter today, I would do it differently, and use an elliptic filter instead of a Butterworth design:

% DEFINE FILTER PARAMETERS 
Fs  = 2E+6;
Fn  = Fs/2;
Fc1 =  630;
Fc2 =  640;
Fs1 = Fc1*0.8;
Fs2 = Fc2/0.8;
Rp  =    1;
Rs  =   10;

% DESIGN FILTER
[n,Wp]  = ellipord([Fc1 Fc2]/Fn, [Fs1 Fs2]/Fn, Rp, Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp);
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^14, Fs) 

The design procedure apparently appears to be important in producing the correct filter result, since even the Butterworth design is significantly impproved:

% DESIGN FILTER
[n,Wn]  = buttord([Fc1 Fc2]/Fn, [Fs1 Fs2]/Fn, Rp, Rs);
[z,p,k] = butter(n,Wn);
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^14, Fs)

I have built a number of hardware filters, however they were all op-amp active filters with discrete components (and none recently), not FPGA filters.

See if this approach produces a better result.

Please post back with a Comment reporting the result you get, and how it works in the FPGA implementation.