Butterworth Bandpass filter design

I'm trying to design a butterworth bandpass filter and the following code is what I have so far:

However, when I run this code my signal_filtered matrix does not represent the correct data. I've set it so the passband frequency is 9kHz-10kHz and it seems like all frequencies are being rejected.

FYI, my text file contains frequency and dB data.

Any help would be greatly appreciated :)

Answer 1

[Star Strider]

A better design would be:

[z,p,k] = butter(n, Wp, 'bandpass');
[sos, g] = zp2sos(z,p,k);

or the complete code (that I typed myself because my version of MATLAB will not run images of code):

Fs = 22500;
Fn = Fs/2;
Wp = [9000 10000]/Fn;
Ws = [8500 10500]/Fn;
Rp = 0.5;
Rs = 50;
[n, Wn] = buttord(Wp, Ws, Rp, Rs);
[z,p,k] = butter(n, Wp, 'bandpass');
[sos, g] = zp2sos(z,p,k);

and you can see the Bode plot:

figure
freqz(sos, 2^16, Fs)

so it is obvious that only a small portion of your signal will be passed. That is how the filter is designed, and it appears to work correctly.

.

[Matija Jankovic]

Thanks for your response Star Strider, I really appreciate it!

I run into the same problems when using the code you provided above. I have attached the text file I'm using for the frequency and dB data.

This first image shows my frequency data in Hz in matlab as the first few lines of my code import:

This second image shows my frequency data in matlab after running the code. As you can see, something weird is going on:

Thanks so much!

[Star Strider]

I am not certain what you are doing with thesde data. They are quite obviously frequency-domain data of some sort (they appear to be the Bode plot of some system), so a time-domain filter is absolutely not appropriate.

fidi = fopen('ThreeSourceSetup2_f.txt');
C = textscan(fidi, '%f(%fdB,%f°)', 'CollectOutput',1);
fclose(fidi);

freq = C{:}(:,1);                           % Hz(?)
ampdB = C{:}(:,2);                          % dB
phsdeg = C{:}(:,3);                         % Phase °

figure
subplot(2,1,1)
plot(freq, ampdB)
ylabel('Amplitude (dB)')
subplot(2,1,2)
plot(freq, phsdeg)
ylabel('Phase (°)')
xlabel('Frequency')

The Bode Plot:

Doing a frequency-domain filter is simply a matter of reconstructing the complex version of these data, annding the complex-conjugate to the end of it, and then doing an inverse Fourier transform to get the approximate time-domain signal.

Example —

absamp = db2mag(ampdB);
phsrad = deg2rad(phsdeg);
cplx_half = absamp.*exp(1j*phsrad);
cplx_whole = [cplx_half; conj(cplx_half)];
Fs = 2*max(freq);
Ts = 1/Fs;
time_domain = ifft(cplx_whole);

t = linspace(0, 1, numel(cplx_whole))*Ts;
figure
subplot(2,1,1)
plot(t, real(time_domain))
grid
subplot(2,1,2)
plot(t, imag(time_domain))
grid

I am not certain what that signal is. Ideally, the inverse Fourier transform of a signal should have a vanishingly-small imaginary part, however this obviously has a substantial imaginary part.

In any event, you can now do time-domain filtering on the reconstructed signal (‘time_domain’ in my code). The filter should work. Note that the sampling frequency ‘Fs’ is 40200, assuming that the maximum frequency in the file is the Nyquist frequency.

.

[Matija Jankovic]

This is actually really quite useful for me Star Strider. Thanks for taking the time to do this! I realize I should have probably provided a little more information about the data itself.

Basically, the text file I shared was the freq/dB data from a circuit I simulated on LtSpice. I have also attached the time domain data text file here which I'm assuming you're wanting to use?

The aim: In this circuit I have 3 AC sources (10kHz, 11kHz and 12kHz). The data in this text file is the circuit response to all three frequencies combined. I am trying to filter out the 11k and 12kHz frequencies responses and only be left with the 10kHz response. Maybe with this new text file I can do this? Thanks again!

[Star Strider]

My pleasure.

If all the sampling intervals are equal, the sampling interval is 3.697882035164080e-007 and the sampling frequency is 2.704250677795428e+06. Use those with the filter deisgn code and you should have no problems. Remember, the second column is the signal, so filter it. Plot the filtered and unfiltered signals against the first column, that ppears to be the sampling times. You have already written most of that (except for the plot calls, assuming that you want to plot the original signal and filtered output), so there is nothing for me to add at this point.

If my Answer helped you solve your problem, please Accept it!

.