From shaehnic@sdcc13 (Steve - Happy Hacker!) Thu Apr 2 19:36:52 1992
From: shaehnic@sdcc13.ucsd.edu (Steve - Happy Hacker!)
Newsgroups: alt.sources
Subject: Simple FFT spectrum analysis package. v1.5
Message-ID: <31248@sdcc12.ucsd.edu>
Date: 2 Apr 92 07:31:41 GMT
Reply-To: shaehnic@ucsd.edu
Organization: HyperHack - San Diego, CA.
Lines: 1470
Nntp-Posting-Host: sdcc13.ucsd.edu
This is a simple but efficient FFT package.
The FFT routines are very modular C that can be used in other
programs, and I have tried to keep things simple and explanatory.
Included is a small time-varying sliding-window audio analysis program
that generates spectrum output suitable for XGraph, Gnuplot, or
further processing.
Enjoy!
-Steve
Steve Haehnichen
E-Mail: shaehnic@ucsd.edu
----------8<----------8<----- CUT HERE -----8<----------8<----------
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# COPYING
# Makefile
# README
# fft.c
# fft.doc
# fft.h
# plot3d.gp
# swindow.c
# swindow.doc
# This archive created: Wed Apr 1 23:23:23 1992
export PATH; PATH=/bin:$PATH
if test -f 'COPYING'
then
echo shar: will not over-write existing file "'COPYING'"
else
cat << \SHAR_EOF > 'COPYING'
GNU GENERAL PUBLIC LICENSE
Version 1, February 1989
Copyright (C) 1989 Free Software Foundation, Inc.
675 Mass Ave, Cambridge, MA 02139, USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
Preamble
The license agreements of most software companies try to keep users
at the mercy of those companies. By contrast, our General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users. The
General Public License applies to the Free Software Foundation's
software and to any other program whose authors commit to using it.
You can use it for your programs, too.
When we speak of free software, we are referring to freedom, not
price. Specifically, the General Public License is designed to make
sure that you have the freedom to give away or sell copies of free
software, that you receive source code or can get it if you want it,
that you can change the software or use pieces of it in new free
programs; and that you know you can do these things.
To protect your rights, we need to make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.
For example, if you distribute copies of a such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have. You must make sure that they, too, receive or can get the
source code. And you must tell them their rights.
We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.
Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software. If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.
The precise terms and conditions for copying, distribution and
modification follow.
GNU GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License Agreement applies to any program or other work which
contains a notice placed by the copyright holder saying it may be
distributed under the terms of this General Public License. The
"Program", below, refers to any such program or work, and a "work based
on the Program" means either the Program or any work containing the
Program or a portion of it, either verbatim or with modifications. Each
licensee is addressed as "you".
1. You may copy and distribute verbatim copies of the Program's source
code as you receive it, in any medium, provided that you conspicuously and
appropriately publish on each copy an appropriate copyright notice and
disclaimer of warranty; keep intact all the notices that refer to this
General Public License and to the absence of any warranty; and give any
other recipients of the Program a copy of this General Public License
along with the Program. You may charge a fee for the physical act of
transferring a copy.
2. You may modify your copy or copies of the Program or any portion of
it, and copy and distribute such modifications under the terms of Paragraph
1 above, provided that you also do the following:
a) cause the modified files to carry prominent notices stating that
you changed the files and the date of any change; and
b) cause the whole of any work that you distribute or publish, that
in whole or in part contains the Program or any part thereof, either
with or without modifications, to be licensed at no charge to all
third parties under the terms of this General Public License (except
that you may choose to grant warranty protection to some or all
third parties, at your option).
c) If the modified program normally reads commands interactively when
run, you must cause it, when started running for such interactive use
in the simplest and most usual way, to print or display an
announcement including an appropriate copyright notice and a notice
that there is no warranty (or else, saying that you provide a
warranty) and that users may redistribute the program under these
conditions, and telling the user how to view a copy of this General
Public License.
d) You may charge a fee for the physical act of transferring a
copy, and you may at your option offer warranty protection in
exchange for a fee.
Mere aggregation of another independent work with the Program (or its
derivative) on a volume of a storage or distribution medium does not bring
the other work under the scope of these terms.
3. You may copy and distribute the Program (or a portion or derivative of
it, under Paragraph 2) in object code or executable form under the terms of
Paragraphs 1 and 2 above provided that you also do one of the following:
a) accompany it with the complete corresponding machine-readable
source code, which must be distributed under the terms of
Paragraphs 1 and 2 above; or,
b) accompany it with a written offer, valid for at least three
years, to give any third party free (except for a nominal charge
for the cost of distribution) a complete machine-readable copy of the
corresponding source code, to be distributed under the terms of
Paragraphs 1 and 2 above; or,
c) accompany it with the information you received as to where the
corresponding source code may be obtained. (This alternative is
allowed only for noncommercial distribution and only if you
received the program in object code or executable form alone.)
Source code for a work means the preferred form of the work for making
modifications to it. For an executable file, complete source code means
all the source code for all modules it contains; but, as a special
exception, it need not include source code for modules which are standard
libraries that accompany the operating system on which the executable
file runs, or for standard header files or definitions files that
accompany that operating system.
4. You may not copy, modify, sublicense, distribute or transfer the
Program except as expressly provided under this General Public License.
Any attempt otherwise to copy, modify, sublicense, distribute or transfer
the Program is void, and will automatically terminate your rights to use
the Program under this License. However, parties who have received
copies, or rights to use copies, from you under this General Public
License will not have their licenses terminated so long as such parties
remain in full compliance.
5. By copying, distributing or modifying the Program (or any work based
on the Program) you indicate your acceptance of this license to do so,
and all its terms and conditions.
6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the original
licensor to copy, distribute or modify the Program subject to these
terms and conditions. You may not impose any further restrictions on the
recipients' exercise of the rights granted herein.
7. The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time. Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.
Each version is given a distinguishing version number. If the Program
specifies a version number of the license which applies to it and "any
later version", you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation. If the Program does not specify a version number of
the license, you may choose any version ever published by the Free Software
Foundation.
8. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission. For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this. Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.
NO WARRANTY
9. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.
10. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
END OF TERMS AND CONDITIONS
Appendix: How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to humanity, the best way to achieve this is to make it
free software which everyone can redistribute and change under these
terms.
To do so, attach the following notices to the program. It is safest to
attach them to the start of each source file to most effectively convey
the exclusion of warranty; and each file should have at least the
"copyright" line and a pointer to where the full notice is found.
Copyright (C) 19yy
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 1, or (at your option)
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
Gnomovision version 69, Copyright (C) 19xx name of author
Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
The hypothetical commands `show w' and `show c' should show the
appropriate parts of the General Public License. Of course, the
commands you use may be called something other than `show w' and `show
c'; they could even be mouse-clicks or menu items--whatever suits your
program.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary. Here a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the
program `Gnomovision' (a program to direct compilers to make passes
at assemblers) written by James Hacker.
, 1 April 1989
Ty Coon, President of Vice
That's all there is to it!
SHAR_EOF
fi # end of overwriting check
if test -f 'Makefile'
then
echo shar: will not over-write existing file "'Makefile'"
else
cat << \SHAR_EOF > 'Makefile'
#
# Makefile for simple-fft. See swindow.c for copying info.
# $Id: Makefile,v 1.6 1992/03/24 01:03:49 steve Exp $
#
CFLAGS = -Wall -O # Use -O2 if you have gcc v2.0
CC = gcc
LDLIBS = -lm
all: swindow
# Time-varying swindows spectrum analyzer
swindow: swindow.o fft.o getopt.o fft.h
$(CC) $(CFLAGS) -o swindow swindow.o fft.o getopt.o $(LDLIBS)
fft.o: fft.h
clean:
rm -f *.o core swindow
SHAR_EOF
fi # end of overwriting check
if test -f 'README'
then
echo shar: will not over-write existing file "'README'"
else
cat << \SHAR_EOF > 'README'
The Simple FFT Package, version 1.5
Steve Haehnichen
April 2, 1992
This package contains some basic tools for FFT sound spectrum analysis.
See the text files "fft.doc" and "swindow.doc" for detailed info.
All parts of this package are covered by the GNU General Public License.
Please read the file COPYING for use and distribution information.
Makefile Standard makefile for swindow.
Be sure to set -O2 if you have gcc2.0, as it results in
a very significant performance increase.
fft.c Basic functions for general FFT use.
fft.h Header file for any user program using fft.c
fft.doc Some introductory info on using FFTs and fft.c.
swindow.c Spectrum analyzer for audio samples.
swindow.doc Detailed explanation of the swindow program.
I look forward to any comments or criticisms of the package, and hope
to enhance both the program and the documentation in the future.
Enjoy!
-Steve
Steve Haehnichen
E-Mail: shaehnic@ucsd.edu
SHAR_EOF
fi # end of overwriting check
if test -f 'fft.c'
then
echo shar: will not over-write existing file "'fft.c'"
else
cat << \SHAR_EOF > 'fft.c'
/*
* fft.c: Code for generic fft routines.
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]
*
* Copyright (C) 1991 Steve Haehnichen.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 1, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*
* fft.c was derived from US Navy radar code written by Steve Sampson
* in June, 1988 for the MS-DOS platform, and placed in the Public Domain.
* The fft() function was copied directly from fft.c in his FFT26.ZIP.
* Steve Sampson
* Box 45668
* Tinker AFB, OK, 73145
*
* Steve Haehnichen (shaehnic@ucsd.edu) March, 1992.
*
* $Id: fft.c,v 1.8 1992/04/02 06:56:34 steve Exp steve $
*/
#include "fft.h"
# define SINE(x) SineTable[(x)]
# define COSINE(x) CosineTable[(x)]
# define WINDOW(x) WinTable[(x)]
static FLOAT *SineTable, *CosineTable, *WinTable;
/*
* This PERMUTE could be defined as a table lookup (Sampson did this),
* but I found that to be significantly slower than just computing it
* inline. Of course, this is w/ GCC, and not MS-DOS Turbo C.
* I'll leave it as a macro; you can tweak with it.
*/
#define PERMUTE(x, y) reverse((x), (y))
/* Number of samples in one "frame" */
#define SAMPLES (1 << bits)
#define REAL(x) wave[(x)].re
#define IMAG(x) wave[(x)].im
#define ALPHA 0.54 /* ala hanning */
/*
* Bit reverser for unsigned ints
* Reverses 'bits' bits.
*/
inline const unsigned int
reverse (unsigned int val, int bits)
{
unsigned int retn = 0;
while (bits--)
{
retn <<= 1;
retn |= (val & 1);
val >>= 1;
}
return (retn);
}
/*
* Here is the real work-horse.
* It's a generic FFT, so nothing is lost or approximated.
* The samples in wave[] should be in order, and they
* will be decimated when fft() returns.
*/
void fft (complx wave[], int bits)
{
register int loop, loop1, loop2;
unsigned i1; /* going to right shift this */
int i2, i3, i4, y;
FLOAT a1, a2, b1, b2, z1, z2;
i1 = SAMPLES / 2;
i2 = 1;
/* perform the butterfly's */
for (loop = 0; loop < bits; loop++)
{
i3 = 0;
i4 = i1;
for (loop1 = 0; loop1 < i2; loop1++)
{
y = PERMUTE(i3 / (int)i1, bits);
z1 = COSINE(y);
z2 = -SINE(y);
for (loop2 = i3; loop2 < i4; loop2++)
{
a1 = REAL(loop2);
a2 = IMAG(loop2);
b1 = z1 * REAL(loop2+i1) - z2 * IMAG(loop2+i1);
b2 = z2 * REAL(loop2+i1) + z1 * IMAG(loop2+i1);
REAL(loop2) = a1 + b1;
IMAG(loop2) = a2 + b2;
REAL(loop2+i1) = a1 - b1;
IMAG(loop2+i1) = a2 - b2;
}
i3 += (i1 << 1);
i4 += (i1 << 1);
}
i1 >>= 1;
i2 <<= 1;
}
}
/*
* Put the samples back in order after the FFT scrambles them.
* (Decimation-in-frequecy)
* Untested and unused because I haven't done IFFT yet.
*/
void
reorder (complx wave[], int bits)
{
int i, j;
complx temp;
for (i = 0; i < SAMPLES; i++)
{
j = PERMUTE (i, bits);
if (i > j) /* So we only do each pair once */
{
temp = wave[i];
wave[i] = wave[j];
wave[j] = temp;
}
}
}
/*
* Scale sampled values.
* Do this *before* the fft.
*/
void scale (complx wave[], int bits)
{
int i;
for (i = 0; i < SAMPLES; i++) {
wave[i].re /= SAMPLES;
wave[i].im /= SAMPLES;
}
}
/*
* Calculate amplitude of component n in the decimated wave[] array.
*/
FLOAT amp (int n, complx wave[], int bits)
{
n = PERMUTE (n, bits);
return (hypot (REAL(n), IMAG(n)));
}
/*
* Calculate phase of component n in the decimated wave[] array.
*/
FLOAT phase (int n, complx wave[], int bits)
{
n = PERMUTE (n, bits);
if (REAL(n) != 0.0)
return (atan (IMAG(n) / REAL(n)));
else
return (0.0);
}
/*
* Initializer for FFT routines. Currently only sets up tables.
* - Generates scaled lookup tables for sin() and cos()
* - Fills a table for the Hamming window function
*/
void fft_init (int bits)
{
int i;
const FLOAT TWOPIoN = (atan(1.0) * 8.0) / (FLOAT)SAMPLES;
const FLOAT TWOPIoNm1 = (atan(1.0) * 8.0) / (FLOAT)(SAMPLES - 1);
SineTable = malloc (sizeof(FLOAT) * SAMPLES);
CosineTable = malloc (sizeof(FLOAT) * SAMPLES);
WinTable = malloc (sizeof(FLOAT) * SAMPLES);
for (i=0; i < SAMPLES; i++)
{
SineTable[i] = sin((FLOAT) i * TWOPIoN);
CosineTable[i] = cos((FLOAT) i * TWOPIoN);
/*
* Generalized Hamming window function.
* Set ALPHA to 0.54 for a hanning window. (Good idea)
*/
WinTable[i] = ALPHA + ((1.0 - ALPHA)
* cos (TWOPIoNm1 * (i - SAMPLES/2)));
}
}
/*
* Apply some windowing function to the samples.
*/
void window (complx wave[], int bits)
{
int i;
for (i = 0; i < SAMPLES; i++)
{
REAL(i) *= WINDOW(i);
IMAG(i) *= WINDOW(i);
}
}
/*
* Undo the window function to restore full magnitude.
* Only works if there are NO zeros in WinTable! (like hanning's)
*/
void unwindow (complx wave[], int bits)
{
int i;
for (i = 0; i < SAMPLES; i++)
{
REAL(i) /= WINDOW(i);
IMAG(i) /= WINDOW(i);
}
}
SHAR_EOF
fi # end of overwriting check
if test -f 'fft.doc'
then
echo shar: will not over-write existing file "'fft.doc'"
else
cat << \SHAR_EOF > 'fft.doc'
* fft.doc: Documentation on the generic FFT routines.
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]
*
* Copyright (C) 1991 Steve Haehnichen.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 1, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
----------------------------------------------------------------------
FFT.DOC
This document explains the general-purpose FFT
routines supplied by fft.c and fft.h, and some theory.
GENERAL THEORY:
If you know nothing about the general idea of an FFT, then this might
not be the best place to try to learn it. See below for a few
excellent texts on the subject, and then play with the code. The
ideas are not difficult, but they are unusual, and not always
intuitive. If you have a fairly good idea of what an FFT is and how
it works, but don't know how to code or use it, then this code is for
you. I have made an attempt to keep the code more explanatory than
efficient, so it should be fairly readable.
It's also perfectly possible to just use the swindow program and
never know a thing about FFTs. Try experimenting with the parameters
until you get something you like, and enjoy. You probably don't want
to read the rest of this document either.
DATA FORMAT:
Every sample or spectrum band is represented as a complex number, with
a "real" part (called re) and an "imaginary" part (called im). For
real sound samples (the kind you will probably be dealing with), the
imaginary component of every sample will be zero, and the real value
will represent the strength of the signal at that point in time. When
representing the spectrum of the sound, each complex number will
represent the real and imaginary parts of the response in that
frequency band. These can then be used to compute the amplitude,
power, or phase of the component at that frequency.
The entire waveform is an array of such complex numbers. The size of
this array MUST be a power of two for the FFT algorithm to work at
all. This is ensured by swindow.c, but fft() does not check the array
for validity. Good sizes for audio data range anywhere from 256 to
65536.
RESOLUTION TRADEOFFS:
Keep in mind that bigger sample chunks will give you both
higher frequency resolution, and lower temporal resolution. Small
sample chunks will give you many windows, and thus high temporal
resolution, but your frequency bands will be further apart. This
tradeoff always exists.
Explained another way, you can use a large amount of samples to find
the frequency content of a very large sample, resulting in very fine
frequency steps in your spectrum. The full frequency range of the
spectrum will always be from 0 to half of your sampling rate (good ol'
Nyquist and all that), but you will have finer steps between the
frequency bands. Unfortunately, this requires a lot of samples.
Since the FFT only really shows you the AVERAGE spectrum over the
whole window of samples, you will then know (with high accuracy) the
spectrum of this large chunk of time.
If you are sampling long, unchanging sounds, this is definitely the
way to go. Another use would be to find the average spectrum of many
sounds, like maybe a human voice. Even though the sound changes its
spectrum many times, you are often interested in the "overall content"
of the voice, and not its characteristics at many specific points in
time.
The other extreme is to use very small windows (i.e. give the FFT only
128 samples at a time). This will enable you to see many windows of a
sound, and get a different spectrum for each section of time.
Unfortunately, your frequency resolution is impaired proportionally,
so you can't make very specific distinctions within your spectrum.
This approach may be good for rapidly changing sounds (like breaking
glass) where you just want to see how it changes over time, and aren't
as concerned with the specific frequencies.
A compromise between these two extremes is to use larger windows (so
you get high frequency resolution), but SLIDE your window over time
slowly, so every sample gets analyzed more than once. This lets you
see changes that occur in a time smaller than your window. This
technique is used by the swindow program to achieve both both good
frequency resolution and good temporal resolution. One of the command
line options is how much overlap to use. An overlap of 0.5 will slide
the analysis window half way each frame. No overlap means to progress
one window size each time. A negative overlap would mean to put gaps
in the analysis (i.e. skip samples), but that is not currently
supported by swindow. Of course, a very large overlap (say, 90%) will
take much more CPU time than little or no overlap. A value of about
.25 is usually good.
USING FFT() IN YOUR OWN PROGRAMS:
The fft() function is very generic, and does nothing more than the FFT
itself. It uses table lookups for all the sine and cosine
calculations, and bit manipulations rather than arithmetic whenever
possible. I highly recommend compiling it with the GNU C Compiler
v2.0, using the -O2 flag for maximum speed optimization. I have
intentionally not hand-optimized parts that this compiler is
intelligent enough to improve. This results in easier to read code
that suffers no performance disadvantage. On a fast 486, it can run
in real-time for most reasonable sampling rates. This can make
real-time spectrum displays possible without a dedicated DSP
microprocessor.
The fft() function itself takes two arguments: the array of complex
samples (as described above), and the power of two that is the window
size. (i.e. If 'bits' = 8, there should be 256 samples in the
analysis.) The sample array will be modified in place, and there is
no return value.
Note that the sample values may be any valid floating point numbers,
and don't have to be realistic at all. Of course, audio data will
rarely have any non-zero imaginary components.
The spectrum resulting form the FFT will always be in bit-reversed
order. This is just an artifact of how the FFT gains its speed. For
example, if you want spectrum component number 212 of an 8-bit
analysis, you would check array element wave[43]. (212 = 11010100,
reversed = 00101011 = 43) If this sounds a bit odd, it is. The macro
PERMUTE() will give you a bit-reversed integer if you need one, but
functions are provided in fft.c to just give you the magnitude and
phase of any component, which is usually what you want.
Just remember that once you do the fft() it is no longer in order, and
probably shouldn't be made so, if you expect to an IFFT on the data
later.
WINDOWING FUNCTIONS:
Any time you just "chop" a piece of sound data out of a bigger sample,
you get a very abrupt start and stop. When played, this sounds like a
Pop at the beginning and end of a sample, which is sometimes
compensated for by good DACs. Unfortunately, this Pop really exists
in the data, and will wreak havoc on any spectrum analysis.
Without going into too much detail, the FFT sees any given collection
of samples as part of a *continuous* waveform. This means that the
FFT will work best if the samples could be played over and over again
in a loop without changing any spectrum characteristics. That pop
will certainly change some characteristics, and usually add some
drastically deviant results.
If we could always chop up our samples in places where there is no
sound at all, we would be set. Unfortunately, this is not usually
possible. The solution is to "window" our samples in such a way as to
reduce the extreme beginning and end samples to zero.
This could be thought of as starting the tape with the volume turned
low, and then slowly turning it up to normal volume. Before stopping
the tape, you slowly turn the volume down to zero again.
Over the average, this leaves the spectrum characteristics
undisturbed, and guarantees that there will be no abruptness or pops
that was not originally in the sample.
Several windowing functions exist, and all have respective strengths
and weaknesses. Sometimes you may want no windowing at all.
A function called window() is provided in fft.c that will apply
a window function that is stored in a lookup table. The supplied
fft_init() will initialize this table to a Hanning function, which is
mostly a raised and inverted period of the cosine wave.
A nice feature of this function is that it is undoable. You can
restore the samples to their original value by calling unwindow().
For analysis only, this is usually not necessary.
SWINDOW:
The best way to explain is by example, so I point you to the swindow.c
code for a simple example of using the FFT code. Further comments on
this program can be found in swindow.doc or the code itself.
REFERENCES:
Bracewell, Ron. 1965. The Fourier transform and its applications.
New York: McGraw-Hill
Moore, F. Richard. 1990. Elements of Computer Music. New Jersey:
Prentice Hall.
Strawn, John. 1985. The Computer Music and Digital Audio Series,
Volume 1: Digital Audio Signal Processing. Los Altos: W. Kaufmann, Inc.
SHAR_EOF
fi # end of overwriting check
if test -f 'fft.h'
then
echo shar: will not over-write existing file "'fft.h'"
else
cat << \SHAR_EOF > 'fft.h'
/*
* fft.h: Defs for fft routines
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]
*
* Copyright (C) 1991 Steve Haehnichen.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 1, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* $Id: fft.h,v 1.6 1992/04/02 06:56:44 steve Exp steve $
*/
#include
#include
#include
/*
* If you want double-precision floats used everywhere,
* define this to be "double"
*/
#define FLOAT float
/*
* Since every math.h seems to have their own idea of what the
* elements in struct complex are called, we define our own.
*/
struct _complx
{
FLOAT re;
FLOAT im;
};
typedef struct _complx complx;
/*
* Functions provided by fft.c
*/
void fft (complx wavetrum[], int bits);
void ifft (complx wavetrum[], int bits);
void scale (complx wavetrum[], int bits);
void window (complx wavetrum[], int bits);
void unwindow (complx wavetrum[], int bits);
void fft_init (int bits);
FLOAT phase (int n, complx wavetrum[], int bits);
FLOAT amp (int n, complx wavetrum[], int bits);
inline const unsigned int reverse (unsigned int val, int bits);
SHAR_EOF
fi # end of overwriting check
if test -f 'plot3d.gp'
then
echo shar: will not over-write existing file "'plot3d.gp'"
else
cat << \SHAR_EOF > 'plot3d.gp'
#
# Example script for plotting swindow spectrum output using GnuPlot
#
# Uncomment these lines for Postscript output.
#set terminal postscript default mono "Times-Roman" 12
#set output "out.ps"
set title "Time-Varying FFT" 0,-2
set ylabel "Freq. (Hz)" 1,-.5
set zlabel "Amplitude"
set xlabel "Time (sec)"
set ticslevel 0.25
#set size 1.2,1.2
set parametric
#set surface
set nokey
set view 120, 30, 1, 1
splot "out" with lines
pause -1 "Done."
SHAR_EOF
fi # end of overwriting check
if test -f 'swindow.c'
then
echo shar: will not over-write existing file "'swindow.c'"
else
cat << \SHAR_EOF > 'swindow.c'
/*
* swindow.c: Perform a time-varying FFT analysis of real samples.
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]
*
* Copyright (C) 1991 Steve Haehnichen.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 1, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* Steve Haehnichen (shaehnic@ucsd.edu) Mar. 1992
*
* $Id: swindow.c,v 1.8 1992/04/02 06:56:17 steve Exp steve $
*/
#include "fft.h"
#include
#include
#define USAGE \
"window_size -o out_file -i in_file -r srate -l overlap -nxgvqw\n\n"\
"window_size: FFT window size. Required.\n"\
"in_file: Input filename (waveform)\n"\
"out_file: Output filename (spectra)\n"\
"srate: Sampling Rate (for axis labeling only)\n"\
"overlap: Window overlap (as a fraction of window size, 0.0 to 1.0)\n"\
"\n" \
"-n: No output\n"\
"-x: XGraph format\n"\
"-g: Gnuplot format\n"\
"-v: Verbose text output\n"\
"-w: Do not apply a windowing function\n"\
"-q: No stderr output (quiet)\n"
/*
* What format the waveform samples are in. (i.e. float)
*/
#define SAMPLE_TYPE unsigned char
/* #define SAMPLE_TYPE float */
/* #define SAMPLE_TYPE short int */
/*
* Output Formats:
* XGRAPH is suitable for a 2-D overlayed plot from most Unix 'graph'
* programs. I use David Harrison's Xgraph.
* Each time-frame is represented as a separate line on the graph.
*
* GNUPLOT is for Gnuplot, of course. Available from prep.ai.mit.edu, etc..
* This format is best suited for a 3-D spectrum-over-time plot.
*
* VERBOSE prints out all kinds of useless information to help me debug.
*
* NO_OUTPUT is a bit-bucket mostly for benchmarking and such.
*/
enum output_format { XGRAPH, GNUPLOT, VERBOSE, NO_OUTPUT };
/* If you don't select an output format, this is what you get */
#define DEFAULT_FORMAT NO_OUTPUT
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
void main (int argc, char **argv)
{
int win_size; /* FFT window size */
int keep; /* How many samples to carry over */
int R = 0; /* Sampling rate (if known) */
int bits; /* bits^2 = win_size */
int format; /* Output format */
int smooth_flag; /* Apply window function? */
int nyquist; /* Nyquist freq. (win_size / 2) */
complx *wave; /* Wave samples/spectrum */
complx *old_wave; /* Portion of last wave to keep */
SAMPLE_TYPE *reals; /* Buffer for user sample data */
FLOAT ffa; /* Analysis Freq. (R / win_size) */
FLOAT overlap; /* Window overlap. 0.0 -> ~1.0 */
FLOAT timestep; /* Seconds per FFT frame */
FILE *Fin = NULL; /* Input stream */
FILE *Fout = NULL; /* Output stream */
int i, got, block;
FLOAT magnitude;
int opt, error_flag; /* GNU getopt(3) variables */
extern int optind;
extern char *optarg;
/* Default options */
format = DEFAULT_FORMAT;
smooth_flag = 1;
overlap = 0.25;
R = 0;
error_flag = 0;
while ((opt = getopt (argc, argv,
"l:i:I:o:O:r:R:qngwxv")) != EOF)
switch (opt)
{
case 'l': /* window overlap fraction */
case 'L':
overlap = atof (optarg);
if (overlap < 0.0 || overlap > 1.0)
error_flag++;
break;
case 'q': /* no stderr output */
fclose (stderr);
break;
case 'i': /* Input file */
case 'I':
if ((Fin = fopen (optarg, "rb")) == NULL) {
perror ("Unable to open data input file");
exit (1);
}
break;
case 'o': /* Output file */
case 'O':
if ((Fout = fopen (optarg, "wb")) == NULL) {
perror ("Unable to open output file");
exit (1);
}
break;
case 'r': /* Sampling rate */
case 'R':
R = atoi (optarg);
break;
case 'w': /* Toggle window smoothing */
smooth_flag = !smooth_flag;
break;
case 'g': /* Output formats */
format = GNUPLOT;
break;
case 'x':
format = XGRAPH;
break;
case 'v':
format = VERBOSE;
break;
case 'n':
format = NO_OUTPUT;
break;
default: /* unrecognized option */
error_flag = 1;
break;
}
if (error_flag || !argv[optind])
{
fprintf(stderr, "Usage: %s " USAGE, *argv);
exit(1);
}
/*
* Compute bits such that 2^bits = win_size
* Error if it's not exact.
*/
win_size = atoi(argv[optind]);
bits = 0;
for (i = win_size >> 1; i != 0; i >>= 1)
bits++;
if (win_size != (1 << bits))
{
fprintf (stderr, "%s: Window size must be a power of two!\n", *argv);
exit(1);
}
nyquist = win_size / 2;
if (R == 0) /* If we don't know the samp_rate, */
R = win_size; /* use the window size. */
/* Open the stdin/out if no data file was named. */
if (!Fin)
Fin = stdin;
if (!Fout)
Fout = stdout;
/* Print output header, if necessary */
switch (format)
{
case NO_OUTPUT:
break;
case XGRAPH:
fprintf (Fout,
"TitleText: Time-Varying FFT\n"
"YUnitText: Amp.\n"
"XUnitText: Hz\n"
"Nolines: False\n"
"Bargraph: False\n");
break;
case GNUPLOT:
fprintf (Fout,
"# Time-varying FFT\n"
"# \n"
"0 0 0.0\n");
break;
case VERBOSE:
fprintf (Fout,
"Window size: %d \tSampling rate: %d\n"
"Smoothing: %d \tOverlap: %4.2f\n",
win_size, R, smooth_flag, overlap);
break;
}
fft_init (bits); /* Rev up that FFT engine */
wave = malloc (sizeof(complx) * win_size);
reals = malloc (sizeof(SAMPLE_TYPE) * win_size);
keep = MIN (abs ((int)(overlap * (FLOAT)win_size)), win_size-1);
old_wave = malloc (sizeof(complx) * keep);
for (i = 0; i < keep; i++)
old_wave[i].re = 0.0; /* 0.0's to copy first time around */
ffa = (FLOAT)R / win_size; /* Fundamental Freq of Analysis */
/*
* Read in each block of samples, do the fft,
* and print the results, until EOF.
*/
for (block = 0; ;block++)
{
/* recall the tail of the last waveform */
memcpy (wave, old_wave, sizeof(complx) * keep);
got = fread (reals, sizeof(SAMPLE_TYPE), win_size-keep, Fin);
if (got == 0)
break; /* done if no more samples */
/*
* Copy in the new samples,
* pad any unfilled space with zeros,
* and save the tail of this wave for next time.
*/
for (i = 0; i < got; i++)
{
wave[i+keep].re = (FLOAT) (reals[i]);
wave[i+keep].im = 0.0;
}
for (i = got; i < win_size; i++)
{
wave[i+keep].re = 0.0;
wave[i+keep].im = 0.0;
}
memcpy (old_wave, wave+(win_size-keep), sizeof(complx) * keep);
/*
* Scale the wave down to size,
* Apply the windowing function,
* and compute the FFT
*/
scale (wave, bits);
if (smooth_flag)
window (wave, bits);
fprintf (stderr, "\rComputing FFT... %d", block);
fft (wave, bits);
/*
* Print the spectrum data.
*/
if (format != NO_OUTPUT)
{
timestep = (FLOAT) (1.0 - overlap) * win_size / R;
switch (format)
{
case XGRAPH:
fprintf (Fout,"\n\"Time %9.5f\"\n", timestep * block);
break;
case GNUPLOT:
fprintf (Fout, "\n");
break;
case VERBOSE:
fprintf (Fout,"\nTime %6.2f\n", timestep * block);
break;
default:
break;
} /* switch (format) */
/*
* Print the first half of the spectrum (0 -> R/2).
* wave[0] and wave[R/2] are special cases, because
* they don't need scaling to be full-strength.
*
* We're just eliminating DC altogether, since it
* really messes up graphs and is generally undesirable.
* For any audio processing, you will need all the components,
* so we don't change the wave[] array at all.
*
* A destructive, but faster method would be:
* wave[nyquist - 1].re /= 2.0;
* wave[0].re = 0.0;
* wave[1].re = 0.0;
*/
for (i=0; i < nyquist; i++)
{
if (i == 0 || i == 1)
{
#ifndef SHOW_DC_ANYWAY
magnitude = 0.0; /* no DC, see above */
#endif
}
else
{
magnitude = amp (i, wave, bits);
if (i != nyquist - 1)
magnitude *= 2.0;
}
switch (format)
{
case XGRAPH:
fprintf (Fout, "%9.1f %9.6f\n",
(ffa * i), magnitude);
break;
case GNUPLOT:
fprintf (Fout, "%9.3f %9.1f %9.6f\n",
timestep * block, (ffa * i), magnitude);
break;
case VERBOSE:
printf ("band:%3d (%6.2f Hz) real:%8.6f imag:%8.6f"
" amp:%8.6f phase:%8.6f\n",
i, (ffa * i), wave[i].re, wave[i].im,
amp (i,wave,bits), phase (i, wave, bits));
break;
default: /* Do nothing for unknown formats. */
break;
} /* switch (format) */
} /* for each band */
} /* if !NO_OUTPUT */
} /* for each frame */
fclose (Fin);
fclose (Fout);
free (wave);
free (reals);
fprintf (stderr, "\rDone. \n");
exit(0);
}
SHAR_EOF
fi # end of overwriting check
if test -f 'swindow.doc'
then
echo shar: will not over-write existing file "'swindow.doc'"
else
cat << \SHAR_EOF > 'swindow.doc'
* swindow.doc: Documentation on the swindow FFT program.
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]
*
* Copyright (C) 1991 Steve Haehnichen.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 1, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
----------------------------------------------------------------------
SWINDOW.DOC
This document explains "swindow", a time-varying spectrum
analysis program.
OVERVIEW:
Swindow is a very simple program. In the future, I may do more
complex analysis or resynthesis, but swindow will always be a simple
example of how to use the FFT algorithm in C to analyze audio spectra.
Swindow takes a stream of raw binary audio samples as input. The
format is specified at compile time by appropriately changing the
line in swindow.c that looks like this:
#define SAMPLE_TYPE unsigned char
Floats, short ints, or doubles are usually used. The Sound Blaster
generates unsigned chars, so this is what I use.
Swindow will analyze one portion of your sample stream at a time (call
this a frame), and output the ordered spectrum of the sound for each
frame. This spectrum data is usually plotted by another program, or
saved to a file.
USAGE:
The general command usage is:
swindow window_size -o out_file -i in_file -r srate -l overlap -nxgvq
There are a lot of options, but most have suitable defaults. Since
the GNU getopt package is used, arguments may be specified in any
order and spacing is not critical. The 'window_size' argument is
always required, and defines how many samples will be analyzed as one
frame.
window_size
The window_size value MUST always be a power of two! This is demanded
by the FFT algorithm, and swindow will abort if an improper value is
given. Some values to try are 256, 1024, and 4096. Smaller values
will result in many more frames, but each frame will take less time.
Larger frame sizes will take longer per frame, but run fewer FFTs.
See the file "fft.doc" for an explanation of the resolution tradeoffs
involved in choosing a window size.
-i in_file
-o out_file
These specify filenames to use for input or output data. If either of
these arguments are not supplied, the program will use stdin and
stdout, respectively.
-r sampling_rate
The sampling rate does not really effect the analysis values at all,
but makes the spectrum ranges correct. Without a correct sampling
rate, the frequency axis labels on your graph will probably be wrong.
If you do not supply a sampling rate, the windowing size is used.
Simply put, the picture of the spectrum will always look the same; if
you specify the true sampling rate, the Hz values will be correct too.
-l overlap
The overlap factor determines how far forward the analysis window will be
moved at each frame. This is expressed as a floating point value
between 0.0 and 1.0. A value of 0 will cause the window to move
forward window_size samples for each analysis. This would analyze
every sample exactly once. A value of 1.0 would be a very poor
choice, since the window would never move along at all, and the
program will run until you stop it, analyzing the same window.
Good values to use are generally 0.25 to 0.75. Of course, more
overlap means more FFTs, which takes longer.
The default value is 0.25.
-w
This flag turns off the Hanning windowing feature, and leaves the
frames entirely at full amplitude. This is generally not desirable,
as it introduces many extra false frequencies into the spectrum.
See fft.doc for a more detailed explanation of windowing functions.
Windowing defaults to being used.
-q
This flag suppresses all output that would normally go to stderr.
This includes a running status message that tells which FFT is
currently in progress. Mostly for use in shell scripts or the like.
Output Formats:
-n No Output
-x XGraph
-g Gnuplot
-v Verbose
Each of these selects a different format for printing the spectrum
data. All printing formats print some header between each frame
stating the start time of that frame.
-n
No Output is very simple, and very fast. It skips the output routine
altogether. It's mostly for just running time benchmarks on fft() and
supporting functions. If you think you see a way to speed things up,
try it out with this option and compare timings.
-v
The Verbose format was generally intended for debugging and
human scrutiny. It prints everything swindow knows about the data.
This is not pretty, and looks something like this:
band: 5 ( 12.21 Hz) real:0.006271 imag:0.002898 amp:2.282987 phase:-0.289631
-x
The XGraph format can be used for plotting the data with most Unix
'graph' compatible programs. The one I use is XGraph, by David
Harrison . It shouldn't be hard to modify
this option to make any other standard graph program work properly.
Since graph is 2D, each time frame is printed on the same display,
right on top of the others. This is probably most suitable when you
want to see an overview of the spectrum characteristics, and are not
concerned with what-happened-when information.
A line of XGraph data looks like this:
0.0 0.000000
With the frequency on the X-axis, and magnitude on the Y-axis.
I frequently pipe the spectrum data straight into XGraph with
excellent results. Note that XGraph can't handle more than 64 frames
in one graph, so chose your sizes carefully!
-g
The Gnuplot format is primarily intended for three-dimensional
contour plots of the spectrum over time. It is not fast to generate,
but very informative for rapidly changing sound samples, like human
voice.
To generate these three-dimensional graphs, redirect the spectrum data
to a file (possibly using the -o swindow option), and try running
Gnuplot using the provided plot3d.gp as a controlling script. Like this:
$ gnuplot plot3d.gp
If you are on a graphic terminal, a plot should soon pop up, or you
can uncomment the lines in plot3d.gp that ask Gnuplot to generate a
Postscript file. If you have many points in your spectrum, or many
time frames, you will need some patience. Complex plots aren't always
as informative as something more simple, either.
A line of Gnuplot output data looks like this:
0.000 12.2 4.565974
In order, that's the start time of that window, the frequency of that
band, and the magnitude for that band.
(That first "0 0 0.0" at the top of the output is to ensure that
Gnuplot does not draw connecting lines between the different window
contours. There does not appear to be any setting in Gnuplot that
does this automatically, so we make sure that the number of samples in
inconsistent in that dimension. Hey, it works.)
TODO:
Boy, that should get a file in itself, but I don't think swindow
should get much more complicated. It is meant to be easy to read and
understand, and hopefully serves that purpose.
I would like to implement these features, though:
- IFFT and simple resynthesis.
This will happen as soon as I find out what the correct coefficient
tables should be for the IFFT. If you know, I'd appreciate mail.
- Negative overlap factors.
These don't work just because of laziness. It should be perfectly
OK to specify a negative overlap percentage and have swindow
skip samples after every processing. For real-time applications
this would probably be necessary. Next release..
- Binary output.
If any other program is actually going to use the spectrum data,
ASCII output is too gross. I'm not sure if I should make some kind
of standard format for the spectrum data, or just dump the bytes with
a preceding count.
- Phase vocoding?
Hey, it's about half there. I just need a good IFFT function, and
a function to unroll the phase. This would likely be left to another
special-purpose program, of course. Many exist already.
- Better comments and such.
If this is going to be a tutorial kind of program, it can never
have too many comments. I can think of a couple simplifications too.
- Better DOCS.
The Net really needs a good introduction document to FFT theory and
related topics. This package is a far cry from totally enlightening.
- An Xwindows back-end.
I would really like to pipe the spectrum data straight into a simple
but fast Xwindows application that would make some crude real-time
graph of the spectrum. Something like XGraph, but faster and less
fancy. Auto-scaling and hardcopies aren't needed. Ideas?
SHAR_EOF
fi # end of overwriting check
# End of shell archive
exit 0