이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Implementation of Narrow band pass filter ( Butterworth)
조회 수: 28 (최근 30일)
이전 댓글 표시
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.
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.
채택된 답변
Star Strider
2015년 3월 19일
편집: Star Strider
2015년 3월 19일
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])
댓글 수: 31
Bharath
2015년 3월 19일
편집: Bharath
2015년 3월 19일
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/182389/image.jpeg)
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
2015년 3월 19일
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
2015년 3월 19일
편집: Bharath
2015년 3월 19일
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.
- 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/182391/image.jpeg)
- 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
2015년 3월 19일
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
2015년 3월 19일
Thanks again for your answers.
- I tried to experiment with different variables mentioned in your #4 point but that doesn't seem to filter more further.
- 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.
- 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?
- 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
2015년 3월 19일
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
2015년 3월 19일
편집: Bharath
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 19일
편집: Bharath
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 19일
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
2015년 3월 20일
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/182398/image.jpeg)
Filtered Orbit ( Please note the axis scale which is 1E+58.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/182399/image.jpeg)
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
2015년 3월 20일
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
2015년 3월 20일
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
2015년 3월 20일
편집: Star Strider
2015년 3월 20일
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
2015년 3월 25일
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
2015년 3월 25일
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
2015년 3월 25일
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
2015년 3월 25일
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
2015년 5월 5일
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?
![](/matlabcentral/answers/uploaded_files/30374/filtered%20signals%20comparsion.jpg)
Star Strider
2015년 5월 5일
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
2015년 5월 6일
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
2015년 5월 6일
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
2015년 5월 6일
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
2015년 5월 6일
‘...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
2021년 2월 18일
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.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Filter Design에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)