Butterworth filter's gain formula does not agree with R's `signal` package

Question

I'm trying to calculate the Butterworth filter gain. If I use the formula mentioned on Wikipedia: $$ G^2(\omega) = \frac{G_0^2} {1+\left(\frac{j\omega}{j\omega_c}\right)^{2n}} $$

I don't get a matching result from calculating the gain directly from the filter's magnitude using R's signal package.

library(signal)

# Butterworth filter


# Gain formula from wikipedia
Butterworth_gain <- function(freq, cutoff_frequency, n = 1) {
        1/(1+(freq/cutoff_frequency)^(2*n))
    
}

bf <- butter(n = 2, W = .6, type = "low")
bfr <- freqz(bf)

plot(bfr$f, abs(bfr$h)^2, ty ='l')
lines(bfr$f, Butterworth_gain(bfr$f,  pi*.6, n = 2), col = 'red')

Answer 1

The function butter() computes the coefficients of a discrete-time ("digital") Butterworth filter, whereas the gain formula you used is valid for a continuous-time ("analog") Butterworth filter. According to the R documentation you can use butter() to compute an analog filter using the plane argument.

The gain of a discrete-time Butterworth filter is obtained by the bilinear transform, which substitutes the "analog" frequency variable by $\tan(\omega/2)$, where $\omega$ is the frequency in radians, normalized by the sampling frequency. For a discrete-time unit gain lowpass filter of order $n$, the squared gain is given by

$$G^2(\omega)=\frac{1}{1+\left(\frac{\tan(\omega/2)}{\tan(\omega_c/2)}\right)^{2n}},\qquad |\omega|\le\pi\tag{1}$$