Summary of responses to Softball Fourier Question
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg755] Summary of responses to [mg591] Softball Fourier Question
- From: mcclure at mugwump.ucsd.edu (Reality's email account)
- Date: Wed, 12 Apr 95 19:12:34 -0700
A few weeks ago, I wrote a message here asking about doing Fourier
Transforms using Mathematica. I got some very helpful replies.
(Thanks!) Given the variety of answers, I'm just going to
append the replies verbatim, without trying to summarize things.
Most of them are pretty concise already.
First, my question was:
>Date: Wed, 22 Mar 95 16:33:51 PST
>From: gysin (To: mathgroup at christensen.cybernetics.net)
>To: mathgroup at christensen.cybernetics.net
>Subject: Softball Fourier Question
>Reply-To: bkeeley at ucsd.edu
>
>
>Hi folks,
>
>I'm a new subscriber who has just begun using Mathematica, and I have what
>is probably a "softball" question, but I can't seem to figure it out. (Note:
>I have read through everything on FT in Wolfram's _Mathmatica_ and have just
>run through several of the tutorials in Evans & McClellan's "Signal Processing
>Packages". However, given my background (philosophy and neurobiology), you
>should feel free to challenge my grasp of basic mathematical concepts. :-)
>
>OK, here goes: I've got an analog recording of five different sinusoidal
>sources of electrical signals. (I'm recording the electrical output of five
>South American electric fish in an aquarium from a differential electrode.)
>These sources are in the range of 100-400Hz. I ran that recording through an
>A/D converter which sampled the signal with a period of 25usec (but I
>eventually want to compare different sampling rates), so I have a data file
>that looks like this:
>
> time signal
> 0.00 -0.770
> 0.03 -0.768
> 0.05 -0.758
> 0.08 -0.753
> 0.10 -0.748
> 0.12 -0.744
> 0.15 -0.731
> 0.18 -0.739
> 0.20 -0.726
>
>...and so on.
>
>I want to recover the fundamental frequencies of the five sources, which I
>should be able to do by running an FT on the data, pulling out the five
>largest peaks in the power spectrum. So far, I've figured out how to
>run "Fourier" on a list consisting the second column of numbers above.
>However, that gives me a spectrum with an unscaled _x_ axis. It shows
>me that there's a large peak at 30 "units", another smaller peak at 42.5
>"units", etc. Here's my question:
>
>How do I factor in the sampling rate in order that the power spectrum I
>get out is in Hertz?
>
>Since this seems like such an obvious thing to do with FTs, I figure
>I must be missing something pretty basic. On other FT programs that
>I've seen, they have an option for rescaling the time dimension (the
>_x_ axis of the signal), but I don't see how to do this in Mathematica.
>Anybody got any ideas?
>
>Thanks in advance,
>
>Brian
>
>
Here are the answers I got:
>From: "NELSON M. BLACHMAN" <blachman at gtewd.mtv.gtegsc.com>
Below is a Fourier-transform package I worked out about five years ago, which
may help you. Mma version 2 hadn't yet appeared. I presented it at the
Mathematica Conference in Redwood City, California, in about January 1991.
--Nelson Blachman, retired from
GTE Government Systems Corp.
Mountain View, California
(* Since most applications of the Fourier transform basically involve
continuous
arguments--say t and f or omega--for the given function h[t] and for its
transform H[f] or H[omega], the discrete Fourier transform should conveniently
approximate the continuous transform, with an easily selected sampling interval
dt for h[t] and number of samples n so as to keep aliasing and smearing from
obscuring the the spectral features of interest. FourierRad[H,n,dt] provides
a means for determining H[omega] with argument in radians per second, for
example, and FourierHz[H,n,dt] yields H[f] with argument expressed, for
example,
in hertz. These can be used to show both the spectrum of any given h[t] and
the
adverse effects of too small a sampling range n dt or too large a sampling
interval dt without any need to alter h[t].
LPR[H,((fmin,fmax),Hmin,Hmax)] and LPI[H,((fmin,fmax),Hmin,Hmax)] display solid
and dashed graphs, respectively, of the real and imaginary parts of H[omega] or
H[f]. PlotRange options are provided, as too small an Hmax, for example, may
cause important features to be lost off the top of the graph. The LPR and LPI
plots can be superposed by means of Show[{%, %%}]. After processing in the
fre-
quency domain, one can return to the time domain by using InverseFourierRad[H]
or InverseFourierHz[H], and the results can be plotted by means of LPR and LPI.
*)
FourierRad::usage="FourierRad[h, n, dt] returns n ordered pairs {omega,
H[omega]}, the values of omega being separated by 2 Pi/(n dt), and 0 being a
middle one. In the computation of H[ ], h[t] is evaluated at n times separated
by dt, with 0 again a middle one. Here H[omega] approximates the integral from
t = -n dt/2 to t = n dt/2 of (2 Pi)^-1 H[t] E^(-I omega t)."
FourierRad[h_,n_,dt_] := Transpose[ { Table[ N[ 2Pi (i - 1 - Floor[n/2])
/( n dt )], {i,n}], RotateRight[ N[( Sqrt[n] dt/(2 Pi))
InverseFourier[ RotateLeft[ Table[ N[ h[(i - 1 - Floor[ n/2 ]) dt] ] ,
{i,n}], Floor[n/2] ] ] ], Floor[n/2] ] } ]
FourierHz::usage="FourierHz[h, n, dt] returns n ordered pairs {f, H[f]}, the
values of f being separated by 1/(n dt), and 0 being a middle one. In the
computation of H[ ], h[t] is evaluated at n times separated by dt, with 0
again a middle one. Here H[f] approximates the integral from t = -n dt/2 to
t = n dt/2 of H[t] E^(-2 Pi I f t)."
FourierHz[h_,n_,dt_] := Transpose[ { Table[ N[ (i - 1 - Floor[ n/2 ])/( n dt
)],
{i,n}], RotateRight[ N[( Sqrt[n] dt) InverseFourier[ RotateLeft[ Table[ N[
h[(i - 1 - Floor[n/2]) dt] ] , {i,n}], Floor[n/2] ] ] ], Floor[n/2] ] } ]
LPR::usage="LPR[H, ((fmin, fmax), Hmin, Hmax)] plots as a solid line the real
part of the output H of FourierRad or fourierHz or its Inverse to the extent
that it lies within the ranges of abscissas fmin, fmax and ordinates Hmin,
Hmax that may optionally be specified--without braces."
LPR[H_, r___] := ListPlot[ Re[H], (PlotJoined -> True,
If[ Length[{r}] == 0, PlotJoined-> True,
PlotRange -> If[ Length[{r}] == 2, {r}, {{#1,#2},{#3,#4}}& [r] ] ] ) ]
LPI::usage="LPI[H, ((fmin, fmax), Hmin, Hmax)] plots as a dashed line the
imaginary part of the output H of FourierRad or FourierHz or its Inverse to
the extent that it lies within the ranges of abscissas fmin, fmax and
ordinates Hmin, Hmax that may optionally be specified--without braces."
LPI[H_, r___] := ListPlot[ Map[{#[[1]],Im[#[[2]]]}&, H], (PlotJoined -> True,
PlotStyle -> Dashing[{0.03, 0.015}], If[ Length[{r}] == 0, PlotJoined ->
True,
PlotRange -> If[ Length[{r}] == 2, {r}, {{#1,#2},{#3,#4}}& [r] ] ] ) ]
InverseFourierHz::usage="InverseFourierHz[H_] returns as many ordered pairs
{t, h[t]} as there are in H, the values of t being those used in FourierHz to
obtain H. Here h[t] approximates the integral from -1/(2 dt) to 1/(2 dt) of
h[f] E^(2 Pi I f t)."
InverseFourierHz[H_] := Block[ {n, dt, HH}, HH = Transpose[ H ]; n = Length[
HH[[1]] ]; dt = - Floor[n/2]/(n H[[1,1]]); Transpose[ { N[ n dt^2 ] HH[[1]],
RotateRight[ N[ Fourier[ RotateLeft[ N[ 1/(dt Sqrt[n]) ] HH[[2]],
Floor[ n/2] ] ] ], Floor[ n/2 ] ] } ] ]
InverseFourierRad::usage="InverseFourierRad[H] returns as many pairs {t, h[t]}
as there are pairs in H, the values of t being those used in FourierRad to
obtain H. Here h[t] approximates the integral from -Pi/dt to Pi/dt of
h[omega] E^(I omega t)."
InverseFourierRad[H_] := Block[ {n, dt, HH}, HH = Transpose[ H ]; n = Length[
HH[[1]] ]; dt = - 2 Pi Floor[n/2]/(n H[[1,1]]); Transpose[ { N[ n dt^2 /( 2
Pi ) ] HH[[1]], RotateRight[ N[ Fourier[ RotateLeft[ N[2 Pi /(dt Sqrt[n]) ]
HH[[2]], Floor[ n/2] ] ] ], Floor[ n/2 ] ] } ] ]
==========================================================
From: Paul E Howland <PEHOWLAND at taz.dra.hmg.gb>
On Wed, 22 Mar 1995, Brian Keeley wrote:
> OK, here goes: I've got an analog recording of five different sinusoidal
> sources of electrical signals. (I'm recording the electrical output of five
> South American electric fish in an aquarium from a differential electrode.)
>
> These sources are in the range of 100-400Hz. I ran that recording through an
> A/D converter which sampled the signal with a period of 25usec (but I
> eventually want to compare different sampling rates), so I have a data file
> that looks like this:
According to Nyquist's sampling theory the minimum sample rate you can
use is twice the bandwidth of the signal. For simplicity, if you assume
you signal runs from 0-400Hz, then the lowest sample frequency you can
use is 800Hz.
> I want to recover the fundamental frequencies of the five sources, which I
> should be able to do by running an FT on the data, pulling out the five
> largest peaks in the power spectrum. So far, I've figured out how to
> run "Fourier" on a list consisting the second column of numbers above.
> However, that gives me a spectrum with an unscaled _x_ axis. It shows
> me that there's a large peak at 30 "units", another smaller peak at 42.5
> "units", etc. Here's my question:
>
> How do I factor in the sampling rate in order that the power spectrum I
> get out is in Hertz?
With the Fourier Transform the frequency resolution is given by 1/NT,
where T is your sample period (in seconds) and N is the number of samples
in the FT. Thus with a sample period of 0.1 seconds and 100 samples,
each frequency bin will be spaced 1/(0.1*100) = 1/10 Hz apart. So that
determines the resolution (or spacing) of the bins. The minimum
frequency given by the FT will be -Fs/2 and the maximum frequency will be
Fs/2, where Fs is the sample frequency (1/T). I can't remember how the
Mathematica FT gives its answer, but very often the results from a fast
FT algorithm are presented as frequencies from 0-Fs/2 followed by those
from -Fs/2 to 0 (rather than -Fs/2 to 0 to Fs/2, which is what you might
expect).
However, given your problem, I wouldn't use a Fourier transform. As you
know you are looking for five sinusoidal tones, you could do the whole
thing more accurately by using a parametric spectral estimator. The
commonly used ones are known as Autoregressive (AR) and Moving Average
(MA) spectral estimators. These model the spectrum by an IIR or FIR
filter (respectively) fed by white noise. If you know you are looking
for five tones, then you would model your spectrum with a five pole
digital filter. An AR model is best at modelling spectral peaks, whereas
an MA model is best at modelling spectral nulls. The advantage of using
an AR or MA estimator is that you can achieve much better frequency
resolution than with a FT. Almost any book on digital signal processing
should point you in the right direction. They are pretty easy to code up
in Mathematica.
Hope this helps - if you have any more questions, please ask.
Paul E Howland
Long Range Ground Radar Systems Section tel. +44 (0)1684 895767
CSS2 Division, Room BY209 fax. +44 (0)1684 896315
Defence Research Agency email: PEHOWLAND at DRA.HMG.GB
Malvern, Worcs, WR14 3PS, UK.
==========================================================
From: "loeffler" <loeffler at arlut.utexas.edu>
Here's a short answer:
If the Fourier transform's "bins" are indexed 0 to N-1 then the
corresponding frequecnies are: sampling_freq*(bin_index/N).
If the Fourier transform's "bins" are indexed 1 to N then the
corresponding frequencies are: sampling_freq*((bin_index-1)/N).
Lastly, since your data is (probably) real valued, then only the
first N/2 "bins" represent valid frequencies.
Good Luck
Charles Loeffler
Applied Research Labs
University of Texas at Austin
==========================================================
From: "Rich Klopp" <rich_klopp at qm.sri.com>
Try using the Mac program Igor. It is ideal for waveform operations like
those you propose. It handles sampling rate units automatically. On the
PC/Windows side, you might try DaDISP or VuPoint. These programs are far
better than Mathematica for dealing with experimental data, and especially
waveform data obtained at uniform sampling rates.
==========================================================
From: beau at vision.arc.nasa.gov (Andrew B. Watson)
The basic fact is that the fundamental frequency (the "unit" frequency)
in the output of the Discrete Fourier Transform is 1 cycle/sequence.
Thus if your data record has a total duration of 0.5 second (for
example) then the unit frequency is 1 cycle/0.5 sec = 2 Hz.
Andrew B. Watson
MS 262-2 NASA Ames Research Center
Moffett Field CA 94035-1000
(415) 604-5419 -3323 fax
beau at vision.arc.nasa.gov
http://vision.arc.nasa.gov/
==========================================================
From: "Ed Boss" <boss at wln.com>
Brian writes:
>How do I factor in the sampling rate in order that the power spectrum I
>get out is in Hertz?
The convention for scaling FFT output is to multiply by 1/(2 delta t) where
delta t is your sampling interval. However, there are some pitfalls in
making spectral estimates based on a single FFT of your data. You might
want to look into the Welch method which averages the results of subsets of
the data. Matlab has a routine for implementing this method in its spectrum
command. Also, see page 553 of Digital Signal Processing by Oppenheim and
Schafer. This is a pretty rigorous treatment of signal processing. For an
introduction, you might want to read The Analysis of Time Series by C.
Chatfield. Chapter 7 discusses spectral analysis.
Does anyone know if Mathematica has an implementation of the Welch method,
perhaps in one of the signal processing packages?
Good Luck,
Ed
----------------------------------------------------------------------
Ed Boss .^ *
NOAA/Pacific Marine Environmental Laboratory / /
/ |
/ \
___.---~ `*------ ...