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);

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.