Butterworth Bandpass filter issues

Hi there! I am trying to make a butterworth bandpass filter in matlab and am having some trouble with it. The specifications for my filter are a center frequency of 10kHz with limits of +/- 30 Hz. In addition to this, I'd like to have 100dB of attenuation at 1kHz away from the center frequency, so if anyone could shed some light on how I could do this that would be awesome!

Right now I am trying to test my filter by simulating a 10kHz sinusoidal signal in LTSpice, exporting the time-voltage data to matlab and then resampling the signal so that it is periodic. I have then made and applied a butterworth bandpass filter with passband 9970Hz - 10030 Hz in the hope that I can completely allow the 10kHz signal to pass through. Right now this is not the case and I am getting some cutoff which progressively increases with time. Running the provided matlab script with the provided text file will illustrate this. I have attached the matlab code as well as the exported time-voltage text file.

Any help would be greatly appreciated!

Answer 1

[Star Strider]

That is likely asking too much of a Butterworth design. Use an elliptic filter instead:

Fs = 44100;                                                         % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency
Wp = [1E+3-30 1E+3+30]/Fn;                                          % Passband Frequency (Normalised)
Ws = [1E+3-100 1E+3+100]/Fn;                                        % Stopband Frequency (Normalised)
Rp =   1;                                                           % Passband Ripple
Rs = 100;                                                           % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp,'bandpass');                             % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

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

Use your own sampling frequency.

.

[Matija Jankovic]

Hi Star Strider, thanks for your response.

I tried the elliptic filter, and the transfer function and filtered signal actually turned out worse than the butterworth design.

Butterworth Filtered Signal

Elliptical Filtered Signal

As you can see in these there is clear distortion of the signal during filtering? Do you have any ideas of how to fix this? The 10kHz filter should be fairly easy just to pass through the band right?

Thanks for your help!

[Star Strider]

The filter is doing what you asked it to do:

D = load('10ktest.txt');
t = D(:,1);
s = D(:,2);

figure
plot(t,s)
grid
title('Original Signal')

[rs,rt] = resample(s,t,2.5E+5);                                     % Resample To Constant Sampling Intervals (Fs = 2.5E+5)
Fs = 1/mean(diff(rt));                                              % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency
Wp = [1E+4-30 1E+4+30]/Fn;                                          % Passband Frequency (Normalised)
Ws = [1E+4-500 1E+4+500]/Fn;                                        % Stopband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 100;                                                           % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp,'bandpass');                             % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

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

rs_filt = filtfilt(sos, g, rs);

figure
subplot(2,1,1)
plot(rt, rs)
title('Original Resampled Signal')
grid
subplot(2,1,2)
plot(rt, rs_filt)
title('Filtered Resampled Signal')
grid

Lrs = numel(rt);
FFTrs = fft(rs)/Lrs;
FFTrs_filt = fft(rs_filt)/Lrs;
TrnsFcn = FFTrs_filt ./ FFTrs;
Fv = linspace(0, 1, fix(Lrs/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(TrnsFcn(Iv)))
grid
title('Filter Transfer Function (From Data)')

As I mentioned, that is asking a lot of any filter. You are taking a very small part of the signal, and the filter is returning the appropriate output. (I needed to resample the signal, since ther was a significant variance in the sampling intervals.) This is not the fault of the filter, and instead the fault of the design.

The elliptic filter Bode plot:

.

There was also an error in the earllier passband and stopband frequencies. Those here are correct.

EDIT — (25 Jul 2020 at 1:20)

As an afterthouight, I considered cascading highpass an lowpass filters instead of using the bandpass design:

Wp = (1E+4-30)/Fn;                                                  % Passband Frequency (Normalised)
Ws = (1E+4-100)/Fn;                                                 % Stopband Frequency (Normalised)
Rp =   1;                                                           % Passband Ripple
Rs = 100;                                                           % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                                 % Elliptic Filter Design: Zero-Pole-Gain 
[sos1,g1] = zp2sos(z,p,k);                                          % Second-Order Section For Stability

Wp = (1E+4+30)/Fn;                                                  % Passband Frequency (Normalised)
Ws = (1E+4+100)/Fn;                                                 % Stopband Frequency (Normalised)
Rp =   1;                                                           % Passband Ripple
Rs = 100;                                                           % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp,'low');                                  % Elliptic Filter Design: Zero-Pole-Gain 
[sos2,g2] = zp2sos(z,p,k);                                          % Second-Order Section For Stability

rs1_filt = filtfilt(sos1, g1, rs);
rs2_filt = filtfilt(sos2, g2, rs1_filt);

figure
subplot(2,1,1)
plot(rt, rs)
title('Original Resampled Signal')
grid
subplot(2,1,2)
plot(rt, rs2_filt)
title('Filtered Resampled Signal')
grid

The result is still apparently unsatisfactory.

You are simply asking too much of the filters.

EDIT — (25 Jul 2020 at 3:3)

Thinking about this further, there appears to be no need to filter this signal, since it alkready more than meets the requirements of any bandpass filter:

What am I missing here?

.

[Matija Jankovic]

Thanks for all your responses Star Strider, I really appreciate it.

I also tried the highpass-lowpass cascading but got the same results as you.

All I wanted to do was test the bandpass filter with this 10kHz signal to see if it was working correctly and whether the signal was passing. Eventually I'm going to be passing a more complex signal with 10k,11k and 12k components which I'll need to separate to see their individual responses.

I might have to think of a different method then.

Thanks for your help though!

[Star Strider]

As always, my pleasure!

The bandpass filters should work, regardless. (I still prefer to elliptic filter implementation.) They will likely perform more appropriately if there are other specific frequencies in the signal. It would also be appropriate to have a wider passband than 0.3% of the center frequency. That is likely too narrow for reliable filter performance, regardless of the type of filter.

[Matija Jankovic]

Sounds great, thanks Star Strider. I've increased it so my passband is now wider from 9700 - 10300 (3% of center frequency). Do you have any suggestions of what I should use for the stopband, passband ripple and stopband attenuation?

Also, I forgot to ask earlier but when you set the stopband attenuation to 100 as you did in your code above, does this account for 100dB of attenuation at 1kHz away from the center frequency?

Thanks!

[Star Strider]

As always, my pleasure!

‘Do you have any suggestions of what I should use for the stopband, passband ripple and stopband attenuation?’

Use whatever makes sense in the context of your signal. I usually porefer a passband ripple of 1 dB and a stopband attenuation of 60 dB. This generally creates stable filters. Excessive stopband attenuation creates very steep transition regions and can produce undesirable results.

’Also, I forgot to ask earlier but when you set the stopband attenuation to 100 as you did in your code above, does this account for 100dB of attenuation at 1kHz away from the center frequency?’

That is the maximum stopband attenuation with respecct to the passband maximum. The passband and stopband frequencies determine the shape of the transition region (between passband to stopband frequencies).

.

[Matija Jankovic]

Thanks Star Strider, that makes sense.

Would you happen to have any code that I could use to plot the FFT of the voltage vs time data?

Thanks for your help!

[Star Strider]

As always, my pleasure!

Here is some generic code (untested here, however similar code has worked before):

t = ...;                                                % Time Vector
v = ...;                                                % Voltage Vector
L = numel(t);                                           % Data Length
Fs = 1/mean(diff(t));                                   % Sampling Frequency
Fn = Fs/2;                                              % Nyquist Frequency
FTv = fft(v)/L;                                         % Fourier Transform (Complex)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:numel(Fv);                                       % Index Vector

figure
subplot(2,1,1)
plot(Fv, abs(FTv(Iv))*2)                                % Plot Amplitude
grid
title('Amplitude')
subplot(2,1,2)
plot(Fv, unwrap(angle(FTv(Iv))))                        % Plot Phase
grid
title('Phase (rad/time unit)')
xlabel('Frequency')

This assumes constant differences between adjacent sampling times (regular sampling frequency).

.

[Matija Jankovic]

Awesome thanks for that. The code works well :)

I'm guessing the yaxis in the amplitude plot, in this case, is in units of volts?

Just out of curiosity, is the following line plotting the peak-to-peak amplitude or just the amplitude? The *2 in the second argument is the reason I ask :)

plot(Fv, abs(FTv(Iv))*2) % Plot Amplitude

Thanks again Star Strider, best help I've ever had with coding and explanations!

[Star Strider]

As always, my pleasure!

Yes. The amplitude is in the same units as the original signal.

That depends on the signal. The amplitude is the absolute amplitude subtracted from the mean of the original signal, so half of the peak-to-peak amplitude. For a simple pure sine curve, this is easy to visualise, for a more complex waveform, very much less so.

The Fourier transform is symmetrical and the length of the original time-domain signal vector,with one half being the complex conjugate of the other half. The energy in the original signal is divided approximately evenly between the two halves. So the multiplication by 2 corrects for that, giving the approximate amplitude of the original signal.

.