The Fugue

I often get confused about the scaling of magnitude and power spectra. I think I’ve worked it out so I’m puttin it here for my own benefit, and maybe the benefit of someone on the internet. If I’m wrong, please tell me.

Let’s take a unit-amplitude sinusoid at 100 Hz, and do a Fourier transform. If we look at the amplitude we’d get two components at 100 Hz and -100 Hz, each with amplitude 1/2. Something like this.

Now, if we do the same with a DFT we will get basically the same thing, except we’ll see the effects of discretization of course. But depending on which variation of the DFT you’re using, the magnitude of the components will be either about 1/2 or about N/2. The difference is, of course, the 1/N factor that you included or left out. In the case of the fft function in Octave (and I assume MATLAB), it’s without the 1/N factor.
max(abs(fft(x,1024)))
ans = 443.23

The power spectrum is the magnitude spectrum squared. abs(fft(x)).^2 or (abs(fft(x))/N).^2. Why would you want to square it? Well there are of course many reasons but the nature of audio signals makes most of the ones I know about moot. Furthermore the relative shape of the power spectrum and the magnitude spectrum is the same.

Try this. Take an audio signal (something more interesting than a single sinusoid) and plot the magnitude spectrum in dB, then plot the power spectrum also in dB, e.g.
plot(20*log10(abs(fft(x))))
figure
plot(20*log10(abs(fft(x)).^2))

They have a different scale, but they have the same shape. If you were doing something like peak-picking, it wouldn’t matter which you used if you’re working in dB.

Tensai asked me how I made my graphs in my ringtones post. I’d like to blather on about the graphs and why they’re cool and how you can make them, because that’s the sort of thing I’m good at.

In the olden days, the spectrogram was invented.

Originally greyscale, they are now usually portrayed in color, with “hot” colors meaning higher amplitude and “cool” colors meaning lower amplitude.

While there are myriad ways to view and/or generate spectrograms, the most convenient for me right now is to do it in Octave. If you’re not familiar with Octave, it’s a MATLAB clone. Octave is libre, MATLAB is insanely expensive. Obviously, I use Octave. I have used MATLAB previously (I had a beta copy, until it left beta), and for the most part they are quite comparable. Octave is a bit slower for some things (less optimized) but I’ve seen Octave outperform MATLAB on some specific tasks. The biggest impedance mismatch is in user interface stuff (MATLAB has sophisticated dialog support) and graphing. Octave has most of the essential graphing functionality (it uses gnuplot to render graphs).

So how do we generate a spectrogram in Octave? First we need to read in the WAV file, then we generate the spectrogram.

[x,sr] = wavread('logchirp.wav');
specgram(x,8192,sr);

8192 is the size of the FFT. I find 8192 is a nice compromise between time and frequency resolution (and computation time), but other powers of 2 (especially 1024) are common as well. Here’s the chirp spectrogram:

Notice how the chirp is logarithmic. To our ears, this sounds like a steadily-rising tone. For this reason when dealing with music it’s often better to look at a log-frequency spectrogram. Otherwise all the low frequencies are scrunched together and the relationships between different pitches (and harmonics) aren’t constant. Here’s a log-frequency spectrogram of the same chirp:

This was generated by the Octave code

logfsgram(x,8192,sr); title('logchirp.wav');

Notice the bleed on the low frequencies, this is because we need a longer FFT in order to get more frequency resolution at low frequencies. This is a tradeoff in time resolution though, and processor time. Experiment with different FFT lengths for extra credit.

To use these functions you’ll need to put the specgram and logfsgram “m-files” in your octave search path (current directory or whatever else you specify in your ~/.octaverc).

I struggled with this on and off for a couple of months. Finally I stumbled on the magic needed to get it working.

Dan Ellis has some MATLAB code for log-frequency spectrograms, but in Octave the graph lacks the custom y tick marks. Here’s the original code:

    yt = get(gca,'YTick');
    for i = 1:length(yt)
        ytl{i} = sprintf('%.0f',logffrqs(yt(i)));
    end
    set(gca,'YTickLabel',ytl);

This code gets the existing tick marks and labels them with the log frequency (instead of the frequency bin).
The problem is that Octave doesn’t populate YTick unless it was manually set already.

Here’s the working code:

    % octave doesn't populate YTick with its tick marks, so we have to set
    % our own. Besides, this way we get nice octave intervals.
    set(gca,'YTick',1:12:M);
    yt = get(gca,'YTick');
    for i = 1:length(yt)
        ytl{i} = sprintf('%.0f',logffrqs(yt(i)));
    end
    set(gca,'YTickLabel',ytl);

The Fastest Fourier Transform in the West has some internal parameters for fine tuning exactly how it performs the FFT. For optimal results, you measure the effect of these parameters directly to find the absolutely optimal way to do that transform, but this is overhead that usually takes longer than just doing the transform the “slow” way, so it is only worth doing when you will be doing a whole bunch of similar transforms (e.g. in a spectrogram, filter, etc.).

So FFTW uses different planners, e.g. “patient” and “estimate”, the latter being the default. No matter which planner you use, FFTW gathers “wisdom” about plans and doesn’t have to recalculate them in the same session. Additionally, this wisdom can be saved to disk (in /etc/fftw/wisdom).

FFTW ships with fftw-wisdom, a program to generate wisdom which you can save to disk for the benefit of all FFTW-using programs. You run it something like fftw-wisdom -c -v -o /etc/fftw/wisdom, which takes several hours. I did this, and I found that the generated wisdom refused to be read. Oops. So I have done something almost as good (and which takes much less time).

In my case, I usually use FFTW within Octave. I’m sure it’s in use many other places (video and audio encoders/decoders, etc.) but when I care about the marginal speed gains it’s usually in Octave. Octave gives you access to the wisdom and planners, so you can save the wisdom to disk. Set the planner to “patient” (or “exhaustive” even):

octave:1> fftw('planner','exhaustive');

Then we need to perform an FFT for each plan that we care to optimize. Octave (and MATLAB) only seem to do complex-to-complex transforms, first transforming real data to complex, so all we need to do is generate a vector and perform the FFT for the various sizes we care about. In my case (audio) this does the trick:

octave:2> x = rand(2^15,1);
octave:3> for i=1:15; fft(x,2^i); end
octave:4> fftw('dwisdom')
ans = (fftw-3.1.3 fftw_wisdom
...
)

Copy that answer string and stick it in /etc/fftw/wisdom. Test it by running fftw-wisdom from the command line. If you do transforms that aren’t powers of 2, or are 2-dimensional, etc., the idea is the same. In the future, Octave will read and add to the wisdom already there, so you can build on it over time. If anyone knows a way to automatically save wisdom on exit from Octave (really, a way to automatically do arbitrary code on exit—the code to save the wisdom is simple enough), that would be a neat trick that I’d love to have.

Oh, and if anyone wants to write a patch for Octave to get it to use a real-to-complex transform when applicable that would be approximately a factor of 2 speedup. I might dive in someday, if the need arises.

In MATLAB, you can do this:

x = 2*pi*[0:100]/100;
figure(1); plot(sin(x));
figure(2); plot(cos(x));
% ...
figure(1); title('sin(x)');
figure(2); title('cos(x)');

In Octave on OS X, I was having trouble. There were two problems. First, the version of Octave I had (2.1.73) had a bug where it would not just add the title “sin(x)”, but replot the current plot (the cosine) and give it the title “sin(x)”. Not what I wanted. I upgraded to 2.9.9 and it now works properly provided you solve the other problem.

The other problem is more a misdocumentation. The docs for figure.m read:

Set the current plot window to plot window N. This function
currently requires X11 and a version of gnuplot that supports
multiple frames. If N is not specified, the next available window
number is chosen.

figure is the user-defined function from the file
/usr/local/share/octave/2.9.9/m/plot/figure.m

But if you look at the source it’s obvious that it works with AquaTerm as well,
but only if the environment GNUTERM is set to aqua. So put this in your
~/.bashrc:

export GNUTERM=aqua

On Linux, running a version of Octave without the bug mentioned above would be
sufficient, as you’re obviously running X11. I have gnuplot version 4.0.