This is some code that I wrote for GNU Octave that calculates the lower cutoff, center frequency, and upper cutoff. From what I understand it should be close or identical to code that would run in MATLAB. Basically a notes perceived pitch is directly proportional to the logarithm of it's frequency. The idea is to pick frequencies so that the difference of the logarithms of two adjacent frequencies are constant throughout the spectrum you want to equalize. It's pretty well commented so even if you don't know Octave you should be able to pick it up.
Edit: I went ahead and made an account but can't comment in reply to n00dles because some kind of reputation thing.
Anyway you're close, n00dles. The difference is my code wants the entire frequency range you want to equalize while the nominal frequency of 20Hz and 20kHz are the center frequencies of the lowest band and the highest band. So since each band is 1/3 of an octave wide, that means the -3dB cutoff of each band is 1/6 of an octave from the center frequency. So the low frequency cutoff for the 20Hz band is 2^(-1/6) * 20Hz = 17.818Hz and the upper cutoff frequency is 2^(1/6)*20kHz = 22.449kHz. So in order to reproduce the ISO 31 band frequencies, you should use log_eq_bands(17.818, 22449, 31).
Yes this means a 31 band equalizer is equalizing sounds you can't hear. Putting the center frequencies at the edges of the audible range allows you to have a flat response over the audible range. A filter has -3db at the cutoff frequency so if that cutoff frequency was at 20Hz, you would have that much attenuation by the time you reached the bottom of the audible spectrum.
I'll go ahead and list the frequencies. Note There seems to be some floating point round off error occurring after the 4th significant figure, so I'll give 5. That's still more than what they put on equalizers.
19.999, 25.170, 31.687, 39.893, 50.224, 63.229, 79.603, 100.22, 126.17 158.84, 1.9997, 251.76, 316.95, 399.03, 502.36, 632.46, 796.23, 1002.4, 1262.0, 1588.8, 2000.3, 2518.2, 3170.4, 3991.3, 5024.9, 6326.2, 7964.4, 10027, 12623, 15892, 20008.
As you can see these numbers match up really well when you round them. There's a .04% error at both the 20Hz, and 20kHz bands which I assume we know with absolute precision. So I would feel safe trusting these figures to be accurate within 0.1%. That's 0.1442% of an octave or 1.73% of a half step on the chromatic scale.
I release this code under the WTFPL License. If you want to credit someone call me bollocks.
##This function is used to design an audio equalizer. The user determines the
##frequency range they want the equalizer to handle and how many bands they want
##the equalizer to have. The function then calculates the upper and lower cut-off
##frequency for each band as well as the middle, a.k.a. nominal frequency.
##The lower bound of the frequencies the user wants the equalizer to handle is
##given in the parameter lowerFreq and the upper bound is similarly given in the
##upperFreq parameter. The number of bands desired is given in the parameter
##numBands. These must all be real positive scalars. numBands will be cast as a 32-
##bit integer if it is not already an integer because fractional bands don't make sense
##The output is a matrix with three rows and a number of columns equal to numBands
##bands(1,:) is the list of lower cut off frequencies.
##bands(2,:) is the list of nominal(center) frequencies
##bands(3,:) is the list of upper cut off frequencies
##This function is intended to be used to design equalizers for audio use so the
##equalizer bands are chosen so the center frequencies sound like the difference
##in pitch is constant. This means
##diff(log(bands(1,:))) == diff(log(bands(2,:))) == diff(log(bands(3,:))) == K
##where K is a constant the upper cut-off frequency of one band will equal the
##lower cut-off frequency of the band above it. ex: bands(1,n+1) == bands(3,n)
##as long as n is an integer such that numBands > n > 0
##the cut off frequencies are designed to sound like they have a pitch in the middle
##of the pitches of the nominal frequencies adjacent to them
##ex: (log(bands(2,n)) + log(bands(2, n+1))) / 2 == log(bands(3,n)) == log(bands(1,n))
function bands = log_eq_bands(lowerFreq, upperFreq, numBands)
bands = 0;
##check for valid inputs
if(nargin() != 3)
print_usage ("get_equalizer_bands(start_freq, end_freq, numBands)");
elseif(!(isreal(lowerFreq) && isscalar(lowerFreq)))
print_usage("lowerFreq must be a real positive scalar numeric.");
elseif(!(isreal(upperFreq) && isscalar(upperFreq)))
print_usage("upperFreq must be a real positive scalar numeric.");
elseif(!(isscalar(numBands) && isreal(numBands) && numBands > 0))
print_usage("numBands must be a real positive scalar numeric.");
elseif(lowerFreq >= upperFreq)
usage("upperFreq must be greater than lowerFreq.");
elseif(lowerFreq <= 0)
usage("lowerFreq must be a positive scalar numeric.");
#cast numBands as an integer if it isn't one already
numBands = cast(numBands, "int32");
##cast frequencies as doubles so we know what precision the following calculations
##are going to be
lowerFreq = cast(lowerFreq, "double");
upperFreq = cast(upperFreq, "double");
##create a bunch of frequencies that are spaced logarithmically
##since every band has a lower, center,and upper frequency and shares the upper
##and lower with the adjacent band, then we need 3 frequencies for the first band
##and 2 for each additional band
frequencies = logspace(log10(lowerFreq), log10(upperFreq), (numBands*2 + 1));
##assign every other frequency to a lower cut-off frequency output
bands(1, 1:numBands) = frequencies(1:2:(columns(frequencies) - 2));
##assign every other freuquency to a nominal frequency output after skipping the first
bands(2, 1:numBands) = frequencies(2:2:(columns(frequencies) - 1));
##assign every other frequency to a upper cut-off frequency output after
##skipping the first two
bands(3, 1:numBands) = frequencies(3:2:(columns(frequencies)));