Unable to filter out a square wave signal using Butterworth filter

Hello all,

I am trying to use Butterworth analog filter to filter out a square wave signal.

Here is my code

clear all; close all;

%create stimuli
f_125 = 125e3;
period=1/f_125;
Nr_periods=1;
t=linspace(0,1,32^2)*period*Nr_periods; 
Fs=numel(t) / (t(end)-t(1));
VinSqu=(square(2*pi*(f_125)*t));

%Create filter
Wp = [500e3 2e6]/(Fs/2); %passband in Hz
Ws = [400e3 100e6]/(Fs/2);%stopband in Hz
Rp = 1; % passband ripple
Rs = 20; % stopband ripple
% Butterworth filter order
[n1,wn1] = buttord(Wp,Ws,Rp,Rs,'s')

% Evaluate filter responses
[zb,pb,kb] = butter(n1,wn1,'s');
[bb,ab] = zp2tf(zb,pb,kb);

[hb,wb] = freqs(bb,ab,32^2);
v_filt=filter(bb,ab,VinSqu);
plot(t,v_filt)

The output (v_filt) looks as follows

There is no filtering being done, although the freqs graph looks fine.

What am I doing wrong here. Any help will be greatly appreciated

Answer 1

[Star Strider]

There are some significant errors in the code. To begin with, if you want to use a digital filter, it is necessary to design a digital filter, not a continuous filter. The coefficients are significantly different, so continuous coefficients will not work as desired in a discrete application. I am not certain what you want from the filter, so I also calculated the Fourier transform of the input signal here, since it appears that the filter eliminates most of the actual input signal. To do the actual filtering use the filtfilt funciton, not filter.

%create stimuli
f_125 = 125e3;
period=1/f_125;
Nr_periods=1;
t=linspace(0,1,32^2)*period*Nr_periods; 
Fs=numel(t) / (t(end)-t(1));
VinSqu=(square(2*pi*(f_125)*t));
%Create filter
Wp = [500e3 2e6]/(Fs/2); %passband in Hz
Ws = [400e3 2.5e6]/(Fs/2);%stopband in Hz           % Corrected Upper Stopband
Rp = 1; % passband ripple
Rs = 20; % stopband ripple
% Butterworth filter order
[n1,wn1] = buttord(Wp,Ws,Rp,Rs)                     % Calculate Discrete — Not Continuous — Filter Coefficients
% Evaluate filter responses
[zb,pb,kb] = butter(n1,wn1);
[sos,g] = zp2sos(zb,pb,kb);                         % Second-Order-Sectioni For Stability

figure
freqz(sos,32^2,Fs);

v_filt=filtfilt(sos,g,VinSqu);                      % Use ‘filtfilt’ For Phase-Neutral Filtering

figure
plot(t, VinSqu)
hold on
plot(t,v_filt)
hold off
grid
legend('Unfiltered', 'Filtered', 'Location','NW')

L = numel(t);
Ts = t(2)-t(1);
Fs = 1/Ts;
Fn = Fs/2;
FTVS = fft(VinSqu)/L;                               % Calculate & Plot Fourier Transform Of Input
Fv = linspace(0, 1, fix (L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTVS(Iv))*2)
grid
xlim([0  5E+6])

This at least creates the correct approach. Experiment to get the result you want.

[RAN]

Hi Star Strider,

Thank you for your reply. In order to clariy your questions

I am actually looking to build an analog filter and not a digital filter because it will be built using analog components on circuit level. Hence I used the 's' in butterord and butter and freqs function.

The signal is at between 100 to 150 kHz continuous time square signal and the filter should filter out everything between 500 kHz to 2 MHz. So the filter should not filter out the fundamental and third harmonic (for the signal in the example stated by me)

I am playing around with your code and will keep this post updated

[Star Strider]

To create an analog filter from the coefficients requires circuit synthesis techniques. The text we used is: Chen, Wai-Kai Passive and Active Filters: Theory and Implementations (Wiley, 1986). There are several other textbooks on similar topics by W-K Chen, so others may be newer and more appropriate to what you want to do. I cannot find specific references to MATLAB in the descriptions, however it would not be difficult to implement those calculations in MATLAB.

[RAN]

Hi Star Strider,

In butter function, if the 's' option is enabled then the filter coefficients correspond to analog filter. I reckon 's' means Laplace domain?. Also to evaluate the filter, freqs function is needed instead of freqz. I have attached a snippet of the code.

wn1=[500e3*2*pi 2e6*2*pi];

% Evaluate filter responses
[zb,pb,kb] = butter(n1,wn1,'bandpass','s');
[bb,ab] = zp2tf(zb,pb,kb);
[hb,wb] = freqs(bb,ab,f*2*pi);

%Plot frequency response
semilogx(wb/(2*pi),mag2db(abs(hb)))
xlim([100e3 max(wb/(2*pi))]);ylim([-100 20])

Now instead of the filter function, to filter out data i have used lsim() to analyze the response in time domain.

[ybu, time] = lsim(butter,VinSqu,t);

So far so good.

Now the poblem i face is with zero phase response, because now i would like to filter the data out (VinSqu) with zero phase, but then the designed filter exhibits non zero phase response.

Btw this is the thread which gave me the idea to analayze in time domain

https://www.kvraudio.com/forum/viewtopic.php?t=252211

Could you please suggest if you have any solution to this problem.

[Star Strider]

To use a filter on discrete data, the filter must be discrete. Analog filter coefficients will not work correctly filtering discrete signals.

To do zero-phase filtering in the discrete domain, use the filtfilt function.

To do zero-phase filtering in the analog domain, use a Bessel filter.

[RAN]

In Matlab example code for Bessel filter, there is a phase shift in the passband.

https://nl.mathworks.com/help/signal/ref/besself.html

Why do you suggest to use Bessel filter for analog?

[Star Strider]

That is not a phase shift. It is the normal phase response. Bessel filters are phase neutral, the reason they are used as anti-aliasing filters. Note that it is not possible to convert an analog Bessel filter to a discrete representation.