# Human ear - freq log scale

I am not sure if it's DSP question or physic.

But I wonder is there any math equation (method) that describe that graph? (It could be anything approximately like that).

I am working on audio freq analyser. And for typical pop song I get graph that Low End is very big but Middle and High End are very low. That's the mathematical reality, but I want to transform it to human ear reality. Could anyone help me?

(Of course I know how to make freq log scale on horizontal axis and dB scale on vertical axis, just have problem to modify it to human ear characteristic)

• there is something called "E-weighting" curve from Robert Wannamaker and Stanley Lipshitz, that will describe the bottom curve (0 dB or "Threshold of hearing") in terms of poles and zeros. it's pretty high order, if i recall, about an order of 50. i found this. May 15 '18 at 23:52
• i guess it was "F-weighting". i lose track of the count. i s'pose someone somewhere has a new-improved "G-weighting" or "H-weighting" somewhere. May 16 '18 at 0:26

So Robert Wannamaker did a real bang-up job parameterizing the 0 dB curve and came up with, what he called F-weighting. It's in this AES paper from the 90s:

you can see that this can be turned into an $H(s)$, and with the bilinear transform, into an $H(z)$.

here is the curve on a dB vs. log/linear freq scale. the log frequency plot is the upside-down version of your 0 dB curve above but is more updated from Fletcher-Munson to the ISO standard which cam from i dunno where.

here's a table of poles and zeros:

• Hey Robert, could you give me one more hint? I am not sure how to use it. I mean $W(f)$. Let's say I have calculated in $dB$ magnitude of some frequency, let's say $-20 dB$. So should I multiply it by $W(f)$ like $-20 * W(f)$ or subtract like $-20 - W(f)$. I tried both but I get results that I am not sure if it's OK. I don't want to speak now about exact results, just can't find how to use f-weighting, that's why I am asking about procedure. Thanks in advance. May 21 '18 at 21:22
• are you using MATLAB? how are you evaluating and displaying this $W(f)$? May 21 '18 at 22:31
• I use C++, I have each freq bin magnitude $M(f)$ expressed as something between 0 and 1. And from that I get magnitude in dB as $M_dB(f) = 20 log(M(f))$ where for $M(f)<=0.004$ I set $M_dB(f) = -68 dB$. So now I am not sure how to use f-weighting. $M_dB(f) = 20 log(M(f)) x f_weight(f)$ or $M_dB(f) = 20 log(M(f)) - f_weight(f)$ . Or maybe anything else??? May 22 '18 at 10:38
• My intuition and logic tells me it should be multiplying, but with subtracting I think my graph looks better. So I am not sure. But of course my error could be somewhere else in algorithm, that’s why first I want to be sure how to use f-weighting, then I can consider where I have mistake in implementation of my freq graph analyser. May 22 '18 at 10:42
• are you using the pole-zero data in table 1 or the $W(f)$ equation in the appendix? May 22 '18 at 10:50

I'm not exactly sure what you are trying to do. Determining the perceived loudness of an arbitrary signal is really quite complicated and typically requires a fairly detailed model of human perception including critical bands, spectral masking, temporal masking, loudness contours, etc.

The equal loudness curves themselves are standardized in ISO 226 https://www.iso.org/standard/34222.html. There are a few Matlab versions available: just google "ISO 226 MATLAB" There is also the loudness toolbox that uses much more advanced models: http://genesis-acoustics.com/en/loudness_online-32.html

It should be comment, but There is not enough characters.

for(float i=low_End_index+1.0f; i<=buffNyquist; i++)
{
kkk++;
double wBefore    = dispLogScale * (log10((i-1.0)*logScaleWidth1) - log10(low_End));
double wCurrent   = dispLogScale * (log10((i-0.0)*logScaleWidth1) - log10(low_End));

if(dataSource->outRealMixed[i] <= 0.004)
logaa = 0.0;
else
logaa = (   ( f_weighting(i) * 20.0 * log10(dataSource->outRealMixed[i]))       +48.0 )  /48.0;

tempCurr+=( logaa  );

if(round(wCurrent) != round(wBefore))
{
fftGraph.lineTo(wCurrent, -(tempCurr * zero_dB * 0.75/kkk) + zero_dB);

tempCurr=0.0;
kkk=0.0;
}
}


and f_weighting(i) :

double GraphAnalyser::f_weighting(int freqIndex)
{
double iTof = (double)freqIndex * 22.049 / (double)buffNyquist;
double z1 = pow(iTof, 2.0);
double z2 = pow(pow(0.58, 2.0) + pow(1.03, 2.0) - z1, 2.0) + 4*pow(0.58, 2.0)*z1;
double z3 = pow(pow(3.18, 2.0) + pow(8.75, 2.0) - z1, 2.0) + 4*pow(3.18, 2.0)*z1;
double p1 = pow(0.18, 2.0) + z1;
double p2 = pow(1.63, 2.0) + z1;
double p3 = pow(pow(2.51, 2.0) + pow(3.85, 2.0) - z1, 2.0) + 4*pow(2.51, 2.0)*z1;
double p4 = pow(pow(6.62, 2.0) + pow(14.29,2.0) - z1, 2.0) + 4*pow(6.62, 2.0)*z1;

double dupsko = ( gg * pow(z1, 3.0)*z2*pow(z3, 3.0) / (pow(p1, 3.0) * pow(p2, 2.0) * pow(p3, 4.0))) * ( pow(pow(10.0, 5.0)/p4, 20.0) );

return dupsko;
}

• f_weighting() looks okay except i am not sure your freqIndex is mapping to frequency right. depends on what buffNyquist is. ( f_weighting(i) * 20.0 * log10(dataSource->outRealMixed[i]) doesn't look right. the f_weighting belongs inside the logarithm and i don't know what the 48.0 is about. and it should be 10 log10() May 22 '18 at 11:53
• Yes I am also not sure :) When you go deeper you will see that source is array (vector) outRealMixed. And buffNyquist is just outRealMixed.size() divide by 2. Maybe I should transpose it to log scale, but as you can see outRealMixed gets straight integer “i” as an argument. And it gives me good result on freq axis, just magnitude problem I have. May 22 '18 at 12:01
• I use buffNyquist as an end of loop because I dont want to calculate and display freq above nyquist freq which are in the fact mirror of the freq below nyquist May 22 '18 at 12:03
• Sorry I made little mess. Concern to freq mapping, as I told it’s I think it’s OK due to fact the source is mapping the same way. Concern to -48.0, it’s just my minimum magnitude value (as you remember magnitude is between 0 and 1, so -48 is for mag less than 0.004. I need to add 48 to show zero value on the screen If mag dB is -48. But if mag dB is 0dB than the screen see it as 48 (max value on display), do you understand what I mean? And why do you claim dB should I calculate with 10 log10 instead 20 log10? Are you sure? I always use 20 log10 for volume sliders, so I thought it’s analogical. May 22 '18 at 13:04
• And what do you mean f_weighting belongs to logarithm? So shoul I use i like that: 10 *log10( outRealMixed[i] * f_weighting(i) ). That is right? (Sorry for formating I write on the phone and see no option to format) May 22 '18 at 13:08