I hadn't been looking at it much but I must say I'm a bit baffled to discover that people are using the bilinear transform to build digital parametric EQ's.
It helps to understand what the bilinear transform is for and what it is not. If you have, say, an analogue low-pass or high-pass filter, the bilinear transform will insure that the relative positioning of the poles and zeros produces the same pass-band behaviour. The entire filter response, from DC to infinity, is mapped to the band from DC to fs/2. This is clearly a non-linear mapping.
If you start with a maximally flat butterworth response, the transformed filter will be maximally flat too. The rolloff will be infinitely compressed near fs/2 though. Of course the passband is also frequency-warped, but you don't see that because it's flat... It does become visible if you make a chebychev filter. The transformed filter will have the same ripple amplitude, but the bumps and troughs will be at different frequencies compared to the analogue original.
The only thing pre-warping does is insure that the corner frequency lands where it should. So only the corner frequency is pre-warped!
This makes the bilinear transform quite unsuitable for designing filters that are supposed to modify the entire frequency response (e.g. speaker response correction), because where do you put the corner? Same for an RIAA correction. It has several corners.
There, the direct transform z=log(s/fs) produces markedly better results. The bumps and troughs land at exactly the right frequencies, but as poles and zeros approach fs/2 the amplitude errors grow larger and larger. If you really need to have it as close as you can, I'd suggest a numerical optimisation routine (maybe someone else knows a closed-form alternative but it certainly isn't common knowledge).