Query about Butterworth Filter

Question:

(a) Design a continuous-time Butterworth lowpass filter with 3 poles. Make sure the magnitude response at the useful tone frequency is equal to 0.99. Plot the magnitude response of the designed lowpass filter and check the specification is met.

(b) Convert the continuous-time lowpass filter to a discrete-time lowpass filter. (You can use either impulse invariance or bilinear transformation.) Plot the magnitude response of the designed lowpass filter and check the specification is met.

(c) Transform the discrete-time lowpass filter to a bandpass filter suitable for detecting the tone at the Teager output. Plot the magnitude response of the designed lowpass filter and check if the filter design is done satisfactorily.

Answer: My MATLAB code for part 1 is:

delta1 = 0.01;

delta2 = 0.005;

Rp = -20*log10(1-delta1);

Rs = -20*log10(delta2);

Omegap = 0.2*pi;

Omegas = 0.3*pi;

%find the order and the natural frequency

[N,Omegan] = buttord(Omegap/pi,Omegas/pi,Rp,Rs);

%finding the filter coefficients

[b,a] = butter(3,Omegan);

[H,Omega] = freqz(b,a,2048);

figure(1)

subplot(2,1,1)

plot(Omega/pi,abs(H));

xlabel('Frequency(\Omega/\pi)'), ylabel('|H(e^{j\Omega})|')

ylim([0 1])

xlim([0 1])

subplot(2,1,2)

plot(Omega/pi,20*log(abs(H)));

xlabel('Frequency(\Omega/\pi)'), ylabel('|H(e^{j\Omega})| (dB)')

ylim([-80 5])

xlim([0 1])

figure(2)

subplot(2,1,1)

plot(Omega/pi,abs(H));

xlabel('Frequency(\Omega/\pi)'), ylabel('|H(e^{j\Omega})|')

title('Class example: Butterworth IIR filter passband')

%limit plot to cover passband region

xlim([0 Omegap/pi])

ylim([1-delta1 1+delta1])

subplot(2,1,2)

plot(Omega/pi, abs(H));

xlabel('Frequency(\Omega/\pi)'), ylabel('|H(e^{j\Omega})| (dB)')

title('Class example: Butterworth IIR filter stopband')

%limit plot to cover passband region

xlim([Omegas/pi 1])

ylim([0 delta2])

figure(3)

pzplot(b,a)

axis square

xlabel('Real'); ylabel('Imag');

title('Poles and Zeros of Butterworth IIR filter')

Problem:

1. Does the code do the work? I calculated N and omegan together but in the next line used N=3 as per the question. Is that a problem?
2. I can't plot the poles and zeros in a square plane. Need solution for that.
3. Can't figure out codes for part 2 and 3. Need help

Answer 1

[Pratik]

Hi Sadi,

As per my understanding, you want to design a three-pole continuous-time Butterworth lowpass filter and visualize its zeros and poles. You then plan to convert it into a discrete-time filter, transform it into a bandpass filter, and finally, plot the magnitude responses of all filters.

The code provided works fine. Since, the number of poles required in the filter is three, using ‘N=3’ should indeed be appropriate.

Please follow the below mentioned steps to convert the low pass filter to a discrete-time filter, then to a bandpass filter, and plot the filter responses.

1. To plot the poles and zeros in a square plane, you can utilize the 'zplane' function. Here's the corresponding code snippet for your reference:

        zplane(b,a); % ‘a’ and ‘b’ are filter coefficients

2. The ‘bilinear’ function can be used to convert continuous-time filter to a discrete-time filter. Please refer to the code snippet below:

        % Sampling frequency        
        fs = 1000; % adjust as needed
        
        % Convert to discrete-time filter       
        [bd,ad] = bilinear(b,a,fs);

3. The ‘lp2bp’ function can be used to convert the lowpass filter into a bandpass filter. Please refer to the code snippet below:

        % Desired bandpass frequencies
        W0 = 500; % center frequency, adjust as needed
        BW = 200; % bandwidth, adjust as needed
                
        % Transform lowpass to bandpass        
        [bb,ab] = lp2bp(bd,ad,W0,BW);  

4. The ‘freqz’ function returns the frequency response of the filter in the form of complex numbers. The 'abs' function can be used to obtain the magnitude, and the 'plot' function can be used to depict the magnitude response. Please refer to the code snippet below:

        % Frequency response of bandpass filter      
        [Hb,Omegab] = freqz(bb,ab,2048,fs);    
        plot(Omegab/pi, abs(Hb));

Please refer to the following documentations of ‘zplane’, ‘bilinear’, ‘lp2bp’, ‘freqz’ and ‘plot’ functions for more information:

• ‘zplane’: https://www.mathworks.com/help/dsp/ref/dsp.allpassfilter.zplane.html
• ‘bilinear’: https://www.mathworks.com/help/signal/ref/bilinear.html
• ‘lp2bp’: https://www.mathworks.com/help/signal/ref/lp2bp.html
• ‘freqz’: https://www.mathworks.com/help/signal/ref/freqz.html
• ‘plot’: https://www.mathworks.com/help/matlab/ref/plot.html

Hope this helps!