Power Spectrum

Date: 01/29/98 at 16:32:49
From: David Grasing
Subject: Power Spectrum

I'm currently writing a program that will show the power spectrum of a 
sound using a microphone as an input (as a 3d graph in real time if I 
can manage it). As of right now my program records samples, then 
performs a FFT on the sampled data. When I do a FFT I get a complex 
number and just take the magnitude. Is there a more efficient way of 
doing this, because I'm only interested in the power not the phase?

Date: 02/01/98 at 11:06:30
From: Doctor Randy
Subject: Re: Power Spectrum

Hi David,

Good question! It shows you're really thinking! I really had to think 
this through myself. 

I don't believe you can gain any efficiency by knowing that you need 
only the magnitude and not the phase. Here's why. The FFT is computed 
as a sum, and there is simply no way to sum complex numbers in polar 
form (magnitude/phase). Therefore you must calculate both the real and 
imaginary parts when performing the FFT, and then derive the magnitude
from these parts in the usual manner.

However, there IS a way to gain efficiency in performing the FFT when 
you know that the input to the FFT is real (as opposed to complex). 
You might suspect this when you realize that the transform of a real,
N-point sequence has what is known as Hermitian symmetry, i.e., 
w[n]=w*[-n], where w* denotes the complex conjugate of w. In more 
practical terms, when you look at the power spectrum result on the 
screen, half of it is redundant!

If we know in advance that the input sequence is real (as your 
microphone signal is), then you can exploit this to make one N-point 
FFT simultaneously calculate the results for two N-point input 
sequences. Let x1[n] be the first N-point sequence, and let x2[n] be 
the second N-point sequence. From these we form the N *complex* input 
values to the N-point FFT as 

  x[n] = x1[n] + jx2[n],  n=0, 1, ..., N-1.

Then take the FFT as usual. The resulting transforms of x1 and x2,
which I denote here as X1 and X2, are calculated from the FFT output 
as

  X1[k] = (1/2) * ( X[k] + X*[N-k] )
  X2[k] = (1/2) * ( X[k] - X*[N-k] ), k=0, 1, ..., N-1.

This comes from _Introduction to Digital Signal Processing_, Proakis
and Manolakis, under chapter 9, "Applications of FFT Algorithms."

If any of this is unclear to you, please write back.

-Doctor Randy,  The Math Forum
 Check out our web site!  http://mathforum.org/dr.math/