Apply a bandpass filter in freq domain

Question

Brad 2020-4-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515736-apply-a-bandpass-filter-in-freq-domain

评论： Star Strider 2020-6-16

Hi. I am new to signal processing. Appreciate any help.

I have a signal with too many data points and high sampleing rate (500 MHz). The main goal is discarding the effect of the data outside a certain freq range (10- 100 kHz) for which our instrument is not calibrated. So I think I need to apply a bandpass filter in freq domain. I am not sure if I do that correctly, since after applying this filter, the magnitude of the data outside this range are still not zero and the filter doesn't discard the data outside range and it also changes the values wihtin the range. What is that I am doing wrong?

Here is the code:

fs = 500e6;
f1 = 10e3, f2 = 100e3;
wn = [f1/fs f2/fs];
[b,a] = butter(20,wn,'bandpass');
[H w] = freqz(b,a,length(x));   % x is the main signal
y = fft(x).*(H);                % Applying the filter

Thank you

Brad

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-4-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515736-apply-a-bandpass-filter-in-freq-domain#answer_424298

在 MATLAB Online 中打开

‘So I think I need to apply a bandpass filter in freq domain.’

Don’t. That is extremely difficult, since the Fourier transform is symmetrical and it would be necessary to do that on the ‘positive’ as well as the ‘negative’ (complex-conjugate) parts, and so is generally not worth the effort.

Assuming that ‘x’ is your time-domain signal, do this instead:

fs = 500e6;
fn = fs/2;
f1 = 10e3, f2 = 100e3;
wn = [f1 f2]/fn;
[z,p,k] = butter(20,wn,'bandpass');
[sos,g] = zp2sos(z,p,k);
figure
freqz(sos,2^14,fs);   % x is the main signal
set(subplot(2,1,1), 'XLim',[0 200E3])                   % ‘Zoom’ Plot
set(subplot(2,1,2), 'XLim',[0 200E3])                   % ‘Zoom’ Plot
x_filtered = filtfilt(sos,g,x);

The transfer-function implementation of that filter is demonstrably unstable. The second-order-section implementation is stable and will likely do what you want it to.

I suggest using a different filter design, perhaps an elliptic filter. See ellipord and ellip as well as zp2sos.

14 个评论
显示 12更早的评论隐藏 12更早的评论

Star Strider 2020-4-6

My pleasure.

I need to clarify some things —

1. The ‘transfer function implementation’ refers to using the ‘b’ and ‘a’ vectors. That is not usually the best way to implement a filter.

2. The ‘2^14’ is the length of the Fourier transform used to calculate the transfer function in freqz. Longer transforms create more detailed results.

3. Yes.

4a. Yes. All subsequent calculations should be done on ‘x_filtered’. The filter design functions take a frequency-domain concept and implement it as a time-domain digital filter.

4b. No. It is already in the time domain, so ifft is not necessary to get it into the time domain. You can use fft if you want to see the frequency components of the filtered and unfiltered signals.

5. Filters in general should only attenuate, not amplify. It may be appropriate to go the traditional route and use buttord first with both the passband and stopband frequencies and atttenuations. That will likely resultt in a filter that does not amplify. I always begin with the order estimation functions and providing all the necessary arguments, then use those results as input to the filter design functions, and then to zp2sos for the final implementation, and of course filtfilt to do the actual filtering.

Star Strider 2020-4-7

在 MATLAB Online 中打开

You’re asking the filters to do what is close to an impossible task.

I designed my own (elliptic) filters, and tried cascading highpass and lowpass filters as well as designinmg a bandpass filter. I cannot get either approach to work effectively with my test signal with those passbands.

Fs = 500E6;
Fn = Fs/2;
x = randn(1, 1E+7);                                     % Create ‘x’
Wp = [10E3  100E3]/Fn;                                  % Normalised Passband (Passband = 25 Hz To 75 Hz)
Ws = [ 7E3  150E3]/Fn;                                  % Normalised Stopband (Passband = 20 Hz To 80 Hz)
Rp =  1;                                                % Passband Ripple/Attenuation
Rs = 100;                                               % Stopband Ripple/Attenuation
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                      % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'bandpass');              % Elliptic Filter
[sos,g] = zp2sos(z,p,k);                                % Second-Order-Section For Stability
x_filtered = filtfilt(sos, g, x);
figure
freqz(sos, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 200E3])                   % Zoom
set(subplot(2,1,2), 'XLim',[0 200E3])                   % Zoom
L = numel(x);                                           % Assumes  x  Is A Vector
FT_x = fft(x)/L;
FT_x_filtered = fft(x_filtered)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:numel(Fv);                                       % Index Vector
TF = FT_x_filtered./FT_x;                               % Transfer Function
figure
semilogx(Fv, 20*log10(abs(TF(Iv))))
hold on
yl = ylim;
plot([1 1]*10E3, yl, '-r', 'LineWidth',1)
plot([1 1]*100E3, yl, '-r', 'LineWidth',1)
hold off
grid

I leave this for you to experiment with.

.

Brad 2020-4-10

编辑：Brad 2020-4-10

Thanks Star,

I experimenteted a lot and could get it to work. I ignored the bandpass and only applied a lowpass filter with second-order as you suggested. I accepted your original answer.

I have a few questions though:

1. So can those parameters such as Rs, Rp, Wp be any values? For instance, through reading online it is mentioned that 20<Rs<120 dB. So Rs can never be smaller than 20? or if the filter works for a smaller Rs value, that should be fine? (for me, Rs=20 worked).

2. The value of Ws and its distance from f2 (i.e., the length of transition region) should also be among a certain range? Or whatever values works is fine? I mean theoritically the transtion region should be as small as possible? or it doesn't mater?

3. Is there a rule of thumb for Ws, something like 10% more than f2?

Thanks again so much for your help

Star Strider 2020-4-10

As always, my pleasure! Thank you!

1. The ‘Rs’ and ‘Rp’ are are the stopband and passband attenuations, so ‘Rs’ must always be greater than ‘Rp’ if the filter is to be effective. I do not recall any specific limits on ‘Rs’. (I usually use 50 for ‘Rs’, however there is no specific value for it.) The ‘Wp’ parameter is the passband edge, and defines the filter. (It, in combination with ‘Ws’ define the transition region.)

2. The relation between ‘Ws’ and ‘Wp’ depend on the filter design and where they are in the spectrum, so in a highpass or bandpass filter this is more of a problem if they are close to 0 than if they are higher in the spectrum. (In that instance, I usually choose ‘Ws’ that is half the value of ‘Wp’ for highpass filters or the lower edge of a bandpass filter.) As a general rule, it is best to avoid short/steep transition regions because this will increase filter complexity and can lead to instability. I generally keep the transition regions to about 5% of the passband edges. This depends on the filter, and the results of preliminary experiments with the filter and the particular signal.

3. The 10% criterion works, however I tend to prefer symmetry, so the same transition regions on both sides of a bandpass or stopband filter. It also depends on the filter, so elliptic filters are more likely to be stable with small transition regions than Butterworth designs.

Again, my pleasure!

请先登录，再进行评论。

Answer 2

Brad 2020-6-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515736-apply-a-bandpass-filter-in-freq-domain#answer_451935

在 MATLAB Online 中打开

Hi Star

I applied the filter this way, as you suggested:

x_lp = lowpass(x, 100E3, 500E6);
x_filtered = highpass(x_lp, 10E3, 500E6);

It works. But the elliptic one works as well:

[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                      % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'bandpass');              % Elliptic Filter
[sos,g] = zp2sos(z,p,k);                                % Second-Order-Section For Stability
x_filtered = filtfilt(sos, g, x);

The results are similar but a bit different. My question is which one is correct or more preferbale? Which one should I pick?

Thank you

Brad

3 个评论
显示 1更早的评论隐藏 1更早的评论

Brad 2020-6-16

How exactly to compare the freqz though?

They both do the filtering similarly. By different result, I meant the final outcome of the ifft was a bit different and I am not sure which one is correct.

2. You mean I should use filtfilt even for bandpass? (or low and highpass)? Because I just used those two lines and it worked.

Thank you

Star Strider 2020-6-16

1. Both are ‘correct’. The one that most closely meets your requirements is the one to use. I cannot determine that.

2. Use filtfilt with the digitalfilter object if you request it (second output). Those functions use it internally to do the filtering, so it is not necessary to do anything else. If you are filtering several signals, it is easier to use the digitalfilter object than to repeatedly run those functions. You can use the save function to save it to a .mat file to load when necessary later.

请先登录，再进行评论。