Apply a bandpass filter in freq domain

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

Answer 1

[Star Strider]

‘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.

[Brad]

Thank you Star for your answer.

I am trying to understand what you recommanded.

1. By "transfer-function implementation", you mean multiplying the fft(x) by the transfer function (line 6 in my code) is unstable? and so you suggest your method?

2. Where does 2^14 come from? and shouldn't the last entry to freqz funtion be the length of the x?

3. What we get eventually, x_filtered, is the filtered signal by the defined filter, but it is in time domain, so you mean (and suggest) the filtering should be done in time domain not the freq domain?

4. If yes, then all of the subsequent calculations in my process, should be done on x_filtered? (including the eventual inverse fft to get the data back into time domain?)

5. I tried to plot both x and x_filtered on the same plot, but x_filtered values are more than 3 orders of magnitude larger. Is there a reason?

Thank you

[Star Strider]

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.

[Brad]

Thank you.

1. Regarding 4b., the whole purpose of my analysis is doing some essential processing in freq domain, so fft and ifft have to be done. Therefore, this method might not be good for me, right?

2. That is why I simply wanted to apply a filter in freq domain by a simple multiplciation. Is there anyway to do that?

2. How can I implement this traditional route and use buttord first with both the passband and stopband?

Many thanks

[Star Strider]

My pleasure.

1. I have no idea why you are doing any signal processing in the frequency domain. That is difficult and time-consuming and not likely to produce the result you want because of problems with truncation effects (Gibbs phenomenon), symmetry, and other complications. If you absolutely must do frequency domain filtering, use the fftfilt function.

2. See 1.

3. See the documentation for the buttord function. The documentation explains it much better than I ever could. (For the record, reading the documentation is how I learned how to use it.)

[Brad]

Thanks!

Final and maybe stupid question, if I manually assign all values (of fft(x) in freq domain) before f1 and after f2 to zero and then take ifft, would that be wrong? Seems too simple to be true to me!

[Star Strider]

My pleasure!

It is too simple to be true! It will be necessary to do that on both the original and complex-conjugate parts of the fft result. One way to do it is to use fftshift, then do the zeroing on symmetrical parts of the shifted fft, then do ifftshift to produce the ‘filtered’ version of the fft, then take the ifft of that result.

I again do not recommend this. Use a time-domain filter and be done with it!

[Brad]

Hi Star

Still trying to implement your method, so this the original x and then x-filtered after your filter. As you see, it is hugly amplified. Is this how it should be?

Thank you

[Star Strider]

It looks like the filter design is wrong. This occurs in situations such as this, where the sampling frequency is 500E6 Hz and the cutoff frequencies are 10E3 to 100E3 Hz. Filtering such a small part of the total frequency spectrum, especially at the low end, will cause problems such as the one you are seeing.

The best approach in this situation might be to cascade a highpass filter with a cutoff frequency of 10E3 and a lowpass filter with a cutoff frequency of 100E3. Fllter the original signal with the first filter, and the output of that filter with the second filter. The order is likely not important. (Design and test those filters separately first.)

[Brad]

How to do that? using lowpass and highpass functions or bandpass and/or butter? Any speific parameters or methods etc. should I use?

Thank you

[Star Strider]

I would do:

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

Try that. It should do what you want.

[Brad]

1. I treid it and there is barely any change. I remember doing this originally which eventually lead to my question, becasue this didn't change anything, so I thought of applying it in freq domain. Any other solution?

2. Could you elaborate on this one:

"It will be necessary to do that on both the original and complex-conjugate parts of the fft result. One way to do it is to use fftshift, then do the zeroing on symmetrical parts of the shifted fft, then do ifftshift to produce the ‘filtered’ version of the fft, then take the ifft of that result".

Thanks a bunch

[Star Strider]

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]

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]

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]

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

[Star Strider]

As always, my pleasure!

To see how they differ, use the second output of the lowpass and bandpass functions to get the ‘digitalfilter’ object and use it with the freqz function (it will instead launch the filter designer function), and use the ‘sos’ of the elliptic filter with the freqz function (that will be exactly that) to see the filter Bode plots. Then, depending on your criteria (speed, efficiency, passband characteristics, etc.) choose the filter that best suits your needs.

Use all the filters with the filtfilt function to do the actual filtering.

[Brad]

1. 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]

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.