double check implementation of Fc correction in butterworth filter

[alcantarar]

Opening this because I closed mine. See https://github.com/alcantarar/dryft/issues/22 for alternative implementation of Winter's cutoff frequency correction and literature sources. I believe correction should be applied to adjusted cutoff frequency Wn, not Fc directly:

C = ( 2 ** ( 1 / n_pass ) - 1 ) ** ( 1 / ( 2 * order ) ) 
Wn = ( np.tan ( np.pi * Fc / Fs ) ) / C  # apply correction factor to adjusted angular cutoff
Fc_corrected = np.arctan ( Wn ) * Fs / np.pi  # Hz
b, a = butter(order, Fc_corrected / Fn)  # produces same coefficients as equations in Winter's book.

[moorepants]

The post that is referred https://www.codeproject.com/Articles/1267916/Multi-pass-Filter-Cutoff-Correction seems to indicate this:

$$ C \tan\left(\frac{\pi}{f_s}f_a\right) = \tan\left(\frac{\pi}{f_s}f_c\right) $$

where $f_c$ is the cutoff frequency, $C$ is the correction factor, $f_a$ is the adjusted cutoff frequency, and $f_s$ is the sampling frequency.

$$ C = (\sqrt{2}-1)^{1/(2n)} $$

[moorepants]

Whereas my code does:

$$ f_a = C f_c $$

$$ C = \frac{1}{(\sqrt{2} - 1)^{\frac{1}{2n}}} $$

[moorepants]

Page 69 of Winter equation 3.8 he defines the cutoff frequency for his butterworth filter formulation as:

$$ \omega_c = \tan{\frac{\pi}{f_s}f_c} $$

[moorepants]

Winter's cutoff frequency definition is for a digital filter and I worked out the correct for an analog filter, but my call to butter() is generating a digital filter. So my correction factor doesn't match for the filter design.

[moorepants]

Analog implementation (incorrect):

    nyquist_frequency = 0.5 * samplerate

    # Since we use filtfilt below, we correct the cutoff frequency to ensure
    # the filter**2 crosses the -3 dB line at the cutoff frequency.
    # |H(w)| = sqrt(1 / (1 + (w / wc)**(2n)))
    # wc : cutoff frequency
    # n : Butterworth filter order
    # |H**2(w)| = 1 / (1 + (w / wc)**(2n))
    # |H**2(wc)| = 1 / (1 + (wc / wa)**(2n)) = 1 / sqrt(2) = -3 dB
    # wa : adjusted cutoff frequency for double filter
    # wa = (np.sqrt(2.0) - 1.0) ** (-1.0 / (2.0 * n))
    correction_factor = (np.sqrt(2.0) - 1.0) ** (-1.0 / (2.0 * order))

    # Wn is the ratio of the corrected cutoff frequency to the Nyquist
    # frequency.
    Wn = correction_factor * cutoff / nyquist_frequency