Design of a digital A-weighting filter with arbitrary sample rate I want to A-weight a time series with arbitrary sample rate. 

An analog A-weighting filter is defined exactly by IEC 61672-1. But there's no definition for a digital filter. One method is to use the bilinear transform (BLT) to convert the analog filter to the digital filter (as done here Applying A-weighting). However this method suffers from extreme warping near nyquist (even when the analog poles/zeros are pre-warped):



Figure 1: A-weighting frequency response comparison where the sample rate is $25600\textrm{ Hz}$. 

Instead I'm thinking of using an algorithm than can design a digital IIR filter with arbitrary frequency response and plugging in the frequency response of the analog A-weighting filter. 


Is this a good approach? 
If so, is there a particular algorithm that would be well suited for this? 


I've looked into MATLAB's yulewalk but I would need a corresponding Python implementation to try out. I've also come across Berchin's FDLS method in a few places, like this question for instance, but all of the links appear to be broken.
It's a common misconception that the approximation of an analog filter by a digital filter must be bad close to Nyquist. This idea might come from the ubiquity of the bilinear transform, for which this is usually indeed the case. Of course, there are certain constraints on the frequency response of discrete-time filters at Nyquist, but they do not necessarily need to result in a bad approximation of an analog filter in that frequency range. The quality of the approximation close to Nyquist depends on several factors, among which are the properties of the frequency response of the analog filter, and the fact whether only the magnitude or also the phase of the analog filter need to be approximated.

I've designed a 6th order IIR filter approximating the frequency response of an analog A-weighting filter, as defined here:

$$H(s)=\frac{k\cdot s^4}{(s+129.4)^2(s+676.7)(s+4636)(s+76655)^2}\tag{1}$$

with $k=7.39705×10^9$.

I chose a sampling frequency of $48$ kHz. The design procedure is a heuristic iterative procedure I came up with some time ago. It's a least squares approximation based on the equation error method, and I might write up all the details some day.

Below is a plot of the design result. Note that both plots go up to Nyquist ($24$ kHz). The left figure shows that one can't see any difference between the logarithmic plots of the analog and the digital frequency responses. The right-hand figure shows the approximation error, defined as the absolute value of the difference of the magnitude responses. The error shows a typical least squares behavior.



You can check out the filter yourself. Here are the coefficients:


b =

   0.169994948147430
   0.280415310498794
  -1.120574766348363
   0.131562559965936
   0.974153561246036
  -0.282740857326553
  -0.152810756202003

a =

   1.00000000000000000
  -2.12979364760736134
   0.42996125885751674
   1.62132698199721426
  -0.96669962900852902
   0.00121015844426781
   0.04400300696788968