Lowpass Butterworth filter design (own poles' calculations)

Question

I'm trying to calculate coefficients for my Butterworth filter of my own. I've found this article that looks really reliable. But I'm trying to implement it in MATLAB and I'm getting bandstop filter. Please find attached matlab script below:

clc, clear all, close all;
%DIGITAL BUTTERWORTH FILTER DESING  
%Butterworth LP with -3dB at fc and N order and G gain
fs = 10; %frequency sampling
fc = 1;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain
wc = 2*pi*fc;

%analog filter design
%poles and zeros   
k = [1:1:N];
%poles searching
%sk = wc * exp(((j*pi)/(2*N))*(2.*k+N-1))
sk = wc * exp(j*(pi/2)) * exp(j*(2.*k+1)*pi/2*N);
z = 0;
G = 1;
[NUM,DEN] = zp2tf(z,sk,G);
%no zeros for LP
Hs = tf(1,DEN)

%filter analysis
figure(1);
subplot(3,1,1)
pzmap(Hs);
subplot(3,1,2)
impulseplot(Hs);
subplot(3,1,3)
stepplot(Hs);    
fvtool(1,DEN);  

I've calculated poles according to equation $10.3.47$ in the article and created transfer function with this poles as a parameters. But fvtool shows bandstop filter instead of lowpass (as described in article). In pzmap I'm getting also pole on right side of $s$-plane so it seems to be unstable. I dont really know how to interpretet this results.

Maks

Answer 1

The poles of a butterworth filter are located in the s-plane at the locations:

$$s_{k} = \sigma_{k} + j \omega_{k}$$

Where:

$$\sigma_{k} = -sin(\theta_{k})$$ $$\omega_{k} = cos(\theta_{k})$$.

And:

$$ \theta_{k} = \frac{(2k+1)\pi}{2N} $$

Where $N$ is the filter order. So, simply create an array of these poles for your given order, and convert them to coefficients of the filter. (See Numerical Recipes in C for a routine for doing this).

When the coefficients are created, perform:

$$ a_{k} \leftarrow \frac{a_{k}}{\omega_{p}^{N-k}} $$

Where $\omega_{p}$ is the passband edge, that was prewarped into the analog domain. Then set the numerator $b_{0}$ to be $b_{0} = \frac{1}{a_{0}}$ for normalization and divide all $a_{k}$ by $a_{0}$, where the $\vec{a}$ are the coefficients of the denominator of course.

What you have now is an analog filter. To get a digital filter, transform the poles and zeros of the denominator and numerator polynomials respectively with the billinear transform.

I don't use MATLAB, so I cant help you with your code, but that is the basic description of the procedure to create a butterworth filter. I have it implemented in C and it works great.

Answer 2

Ok, I've solved this problem myself. I've used equation from this article mentioned in first post. I've just implement it in matlab not enough precise.

clc, clear all, close all;

%ANALOG BUTTERWORTH FILTER DESING
%https://engineering.purdue.edu/~ee538/AnalogFilterDesigns.pdf
%Butterworth LP with -3dB at fc and N order and G gain

wc = 1000;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain 

%analog filter design
%poles and zeros 

%poles searching
for k = 0:N-1
   sk(k+1) = wc*exp(j*(pi/2))*exp(j*(2*k+1)*(pi/(2*N)));
end

[NUM,DEN] = zp2tf([],sk,G);
b0 = DEN(end)*G      %gain normalization  

Hs = tf(b0,DEN)

figure();
freqs(b0,DEN)
figure();
plot(sk./10^2, 'x')
xlim([-10,10]);
title('Poles placement in s-plane')

Many thanks for your help, Maks