New here and I just wanna ask the question of what happens when your pre warping doesnt work on your bilinear transform?

I am assuming my isnt working as part of my research I have read that the higher the frequency is the more non linear it becomes, so I am assuming that my frequency I chose just has to many non linearities that it wouldn't ever get re-aligned to the continuous function critical frequency. It tries but doesnt get there.

Sampling Frequency = 44.410kHz

Nyquist Frequency = 22.205kHz

Critical Frequency = 20kHz

Ratio = Fs/Nyquist frequency = 20/22.205 = 90.1%

I noticed the closer you are critical frequency is to your Nyquist frequency when applying the bilinear transform it gets ugly. And it looks like the pre warping can only do so much.

The only way I see around this if I just changed my Sampling frequency to like 96kHz which would make my Nyquist frequency = 48kHz would be then the ratio = 42% this would ease the constraint of the non-linearities and have the pre warping affect work much better, but this raises the question is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz? Because the way this is going I dont see it.

enter image description here


s = tf('s');

Fs = 96e3;
Ts = 1/Fs;
Fc = 20e3;

GP = -3;
GS = -20;
WP = 2*pi*Fc;
WS = 251327.412;

n = ceil(log10((10^(-GS/10)-1)/(10^(-GP/10)-1))/(2*log10(WS/WP)));

WC_1 = WP/((10^(-GP/10)-1)^(1/(2*n))); 
WC_2 = WS/((10^(-GS/10)-1)^(1/(2*n))); 

LP = 1/((s/WC_1)^4+2.613*(s/WC_1)^3+3.414*(s/WC_1)^2+2.613*(s/WC_1)+1);
LP_2 = 1/((s/WC_2)^4+2.613*(s/WC_2)^3+3.414*(s/WC_2)^2+2.613*(s/WC_1)+1);

LP_44_410kHz = c2d(LP,1/44.410e3,['Method','tustin','PrewarpFrequency',WP]);
LP_96kHz = c2d(LP,Ts,['Method','tustin','PrewarpFrequency',WP]);

options = bodeoptions;
options.FreqUnits = 'Hz';
options.Title.String = 'With Pre-Warping';
  • $\begingroup$ boy, that's a strange sample rate. i have never seen it before. are you sure it's not 44100 Hz? $\endgroup$ Sep 24, 2020 at 18:27
  • $\begingroup$ you're aware of the Audio EQ Cookbook, no? $\endgroup$ Sep 24, 2020 at 18:31
  • $\begingroup$ Bwahah it is a strange sample rate thats for sure, I have not $\endgroup$
    – Pllsz
    Sep 24, 2020 at 18:54

2 Answers 2


Regarding the prewarping, it gets significantly challenging to match the analog shape as we approach the Nyquist boundary since the bilinear transform warps infinity in s to z=-1 in z (Nyquist).

Prewarping implies the filter design is being mapped from s to z which wouldn't be my first choice for filter design unless I specifically needed to "copy the analog", which we should avoid given all the tools and capabilities today operating directly on the math of delay and multiplies rather than the limitation of analog components. The least squares algorithm is often my first go-to choice for standard FIR filter design with results as follows:

First the number of taps needed for an FIR solution is driven by the transition bandwidth (in contrast to the ratio of the cut-off relative to the sampling rate). So here for this application that ratio of concern would be $$\Delta f/(f_s) = f_\Delta = (22.205-20)/44.410 \approx 0.05$$ assuming our requirement is to achieve a rejection target by the Nyquist boundary. (Importantly, this is for the example of it specific to the OP's question but the bigger concern is the analog filter that would be required prior to the ADC-- the sampling ratio drives the requirements on that filter to reject aliasing that would not be resolvable digitally, so there may be motivation to increase the sampling rate to simplify that filter- more on this at the end).

Continuing... given this ratio we can use various "back of the napkin" estimators to determine the approximate number of taps neeeded to get 50 dB rejection (and importantly notice that it IS this ratio that drives that):

fred harris: $N \approx \frac{A}{22 f\Delta} = 50/(22 \times 0.05) = 45$

Kaiser's empirical formula: $N \approx (A-8)/(14.4 f_d) = (50-8)/(14.4 \times 0.05) = 58$

Let's see which is closer to meeting the mark using the least squares algorithm (firls in Octave/Matlab/Python scipy.signal):

Freq Response

Zoom in

As the DSP variant Goldilocks and the three filters, in this case fred harris estimator is too little, Kaiser's is too much and as is typical these are estimators, starting points, from which we iterate to find the final solution which in this case was 49 taps (closest to fred). Note I used python scipy.signal.firls which requires the number of taps to be odd; which is why I implemented the Kaiser estimate with 57 taps instead of the estimated 58, and that 48 taps is likely the closest least squares solution given the targets provided. Also as we see by the estimators, doubling the sampling rate would have the advantage of cutting the number of filter coefficients required in half (or even much more depending on what the actual rejection needed) to achieve the same performance.

Finally again note since it is critical; such a sampling rate and desired passband implies (assuming this is the ADC sampling rate as well) an analog filter that would pass DC to 20 KHz with an anti-alias rejection requirement starting at $2 \times 22.205 - 20 = 24.41$ KHz. As for how this effects analog filter complexity and Q's required which ultimately limit this, this is a ratio of $4.41/22.205 \approx 0.2$ which is not overly challenging in the analog world but in simplest form may introduce group delay distortion that can be compensated for digitally. I am not really sure what the purpose of this filter is; if the analog filter is sufficient then this filter may actually not be doing anything. A typical approach is to relax the requirements on the analog filter to only reject the images and then do the detailed filtering digitally, which in this case would be achieved with a higher sampling rate.

  • $\begingroup$ Appreciate the comment, however you're touching a topic I never even touched yet (FIR). I am more familiar with the IIR and converting it to into the S -> Z domain with the classic tustin method. It looks like creating a new filter with FIR would make it work however I am trying to get this little get off the floor with the knowledge I know. But I am glad its possible with an alternative method $\endgroup$
    – Pllsz
    Sep 24, 2020 at 17:11
  • $\begingroup$ This filter purpose is a low pass filter for audio that has a cut off frequency of 20kHz. $\endgroup$
    – Pllsz
    Sep 24, 2020 at 17:12
  • $\begingroup$ @Pllsz Doing it in IIR has significant challenges including not being linear phase, possible instability, quantization issue that only get more challenging as the order of the filter increases. If your goal is this result the FIR approach is so simple in comparison. Please pay attention too to the points I made about the analog filter for if you are sampling at 44.4 KHz to get a 20 KHz passband, most if not all of your filtering would have already been done in the analog (and perhaps for this reason it isn't done adequately in which case you should go to a higher sampling rate). $\endgroup$ Sep 24, 2020 at 17:17
  • $\begingroup$ and if you can implement it with a FIR, the coefficents are determined simply by typing coeff = firls(49, [0, 2*fp/fs, .99, 1], [1, 1, 0, 0] in python (similar form in octave and Matlab) and your done. As a linear phase filter this requires 25 multipliers since the coefficients are symmetric. Pretty simple so a compelling tool for your toolbox rather than coping the analog. $\endgroup$ Sep 24, 2020 at 17:21
  • 1
    $\begingroup$ Just maybe, alright I guess I can close this thing. Pleasure talking to you. So TL;DR increase the Fs rate or change over to FIR $\endgroup$
    – Pllsz
    Sep 24, 2020 at 18:10

is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz?

Yes it is possible... . I use MIM (Magnitude Invariance Method) for this. Here's plot showing 2nd order LPF responses for fc=20kHz-25kHz (1kHz step) @ 44.kHz sampling:

enter image description here

MIM method works well for even lower sample rate implementations but, there are also issues not sorted out yet.

Here's Matlab implementation of the MIM and PIM (Phase Invariance Method) method.

  • $\begingroup$ How did you know to use that method and not the tustin, and why would one use that method. Ontop of that do you know why the tustin method isnt working? Are my concerns correct? $\endgroup$
    – Pllsz
    Sep 24, 2020 at 6:30
  • $\begingroup$ I see thats an interesting conversation that I have never seen before actually. Not a problem I do appreciate the comment $\endgroup$
    – Pllsz
    Sep 24, 2020 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.