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 / / / | / \ ___.---~ `*------ ...

**Re: Calling functions before loading the package**

**ListSphericalPlot???**

**Precision & Accuracy: SUMMARY**

**ListSphericalPlot???**