Message Boards Message Boards

3
|
16409 Views
|
17 Replies
|
14 Total Likes
View groups...
Share
Share this post:

Construct a function from "Fourier" - data

Posted 9 years ago

Hi all,

I have a question which for sure has been answered a 1000 times before, but I didn’t find a satisfactory and convenient solution, albeit I have browsed the web and the Wolfram- (and other ) sites for hours.

I have a list of experimental data. It is known that there is some periodicity (or even more than one). Now I apply “Fourier” which produces a nice other list of complex numbers. When plotting the absolute values of this list (a spectrum) I even find some “peaks” and I think I translate them correctly to frequencies. Now I try to construct a function which approximates my original data, and fail drastically. So I set out to get insights with a “known” example. But I failed as well. With some hints of Daniel Lichtblau (thanks, Dan) I could construct a procedure which worked for my simple example, but it failed already for a slightly more complicated case ( and for sure for my original data ). See the attached notebook.

I am looking for a fool-proofed (the fool am I) method to construct a good approximation for the original function, or better the values given.

If anybody is interested in the original problem I will of course send all data I have.

Kind regards Hans

P.S. I am using version 7

Attachments:
POSTED BY: Hans Dolhaine
17 Replies

I asked Wen-Shin Lee about this and she suggested that it is a natural for Prony's method. This is an area of expertise for her so I decided to attempt an implementation. I will remark that it is a close cousin to sparse polynomial interpolation methods of Ben-Or and Tiwari, and Zippel, and appears in numerous guises in the literature.

What follows is not necessarily best possible but it should be at least reasonably stable numerically (uses singular values under the hood, in two places). For signals with noise it can be modified to behave, e.g. by using a loose tolerance when computing the initial matrix rank.

I will illustrate with the example that has three frequencies (one of which is zero, that is, there's a constant term). The code and expositions I have looked at use time units of 1 so I have a timescale value to compensate for that. Also I show some intermediate results to give some idea of what is being calculated. But really an understanding of this will require some delving into the literature on this method.

func3freq = 0.3 + 3.5 Cos[2.6 t + .8] + 1.2 Cos[5.3 t - .4];
obstime = 80;
timescale = 1/10;
data3 = Table[func3freq, {t, 0, obstime, timescale}];
len = Length[data3];

Set up initial Prony matrix. We won't actually use it, we just need the rank. That tells us the numer of exponentials we are going to fit.

mat = Most[Partition[data3, Floor[len/2], 1]];
Dimensions[mat]
keep = Length[SingularValueList[mat]]

(* Out[321]= {401, 400} *)

(* Out[322]= 5 *)

Now that we know the rank we'll set up the actual matrix used to determine frequencies.

mat2 = Most[Partition[data3, keep, 1]];
Dimensions[mat2]
rhs = Drop[data3, keep];
Length[rhs]

(* Out[334]= {796, 5} *)

(* Out[336]= 796 *)

Find the Prony polynomial coefficients. Then solve for the complex exponentials as roots of the Prony polynomial.

soln = PseudoInverse[mat2].rhs
roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx]
Map[Norm, roots]
freqs = Log[roots]

(* Out[341]= {1., -4.6583940973, 8.99362652134, -8.99362652134, \
4.6583940973} *)

(* Out[342]= {0.862807070515 - 0.505533341205 I, 
 0.862807070515 + 0.505533341205 I, 0.966389978135 - 0.257080551892 I,
  0.966389978135 + 0.257080551892 I, 1.} *)

(* Out[343]= {1., 1., 1., 1., 1.} *)

(* Out[344]= {-1.02917674383*10^-13 - 0.53 I, -1.02917674383*10^-13 + 
  0.53 I, 3.16698056668*10^-13 - 0.26 I, 
 3.16698056668*10^-13 + 0.26 I, -2.40141240226*10^-13} *)

Now do a least-squares fit to get the amplitudes. Recover the trig polynomials. Massage to put into the expected form.

newmat = Map[roots^# &, Range[0, obstime, timescale]/timescale];
Dimensions[newmat]
coeffs = PseudoInverse[newmat].data3
newf = Chop[TrigExpand[ExpToTrig[coeffs.Exp[freqs*t/timescale]]]]
Expand[TrigFactor[newf]]

(* Out[351]= {801, 5} *)

(* Out[352]= {0.552636596429 + 0.233651005386 I, 
 0.552636596429 - 0.233651005386 I, 1.21923674127 - 1.25537315885 I, 
 1.21923674127 + 1.25537315885 I, 
 0.300000000025 + 1.38460601482*10^-18 I} *)

(* Out[353]= 0.300000000025 + 2.43847348254 Cos[2.6 t] + 
 1.10527319286 Cos[5.3 t] - 2.5107463177 Sin[2.6 t] + 
 0.467302010771 Sin[5.3 t] *)

(* Out[354]= 0.300000000025 + 3.49999999956 Sin[2.37079632674 + 2.6 t] + 
 1.20000000005 Sin[1.17079632681 + 5.3 t] *)

This is an essentially perfect recovery of the original signal.

POSTED BY: Daniel Lichtblau

Hello Dan,

thanks again . I will have a look at all these things.

Kind regards

Hans

POSTED BY: Hans Dolhaine

Hello Henrik,

here are my final remarks to your final remarks.

First of all thank you again for your remarks and the work you spent on the problem.

But : of course there are two and only two frequencies in the problem. It is constructed exactly in this way. We could think of some physical process known to have two periodicities and described by that "Law". The task is to find the parameters of the "Law", namely frequencies , amplitudes and phases. And indeed, in the "spectrum" we see two peaks. Great.

So, when you say there are lots of frequencies these are artefacts of both, observation and / or computation (finite observation interval, finite observation time and so on). If you will that is in fact a technical limitation and not a deep principle of nature,

So, what I learnt here to get what I want is to plot a spectrum, get the number of peaks, extract the frequencies as exactly as possible and do a Fit for amplitudes and phases, and perhaps the frequencies as well.

Concerning the book about Fourier.... eons ago I read Ziessow "On-line Rechner in der Chemie" with some chapters about FFT, and Lighthill, "Fourier-Analysis", the latter not being really practical. But I feel these books are not of much use for my very practical problem.

Thanks a lot again and Kind regards

Hans

POSTED BY: Hans Dolhaine

Another approach is to work with a continuous periodogram and try to isolate frequency peaks from that. I'll demonstrate with the running example involving two frequencies.

ff2 = 0.3 + 3.5 Cos[2.6 t + .8] + 1.2 Cos[5.3 t - .4];
obstime = 80;
data2 = Table[ff2, {t, 0, obstime, .1}];
normdata2 = data2 - Mean[data2];

pgram[w_?NumberQ, vals_, gap_] := 
 Abs[Total[vals*Exp[I*w*gap*Range[0, Length[vals] - 1]]]^2]/
  Length[vals]

gap = obstime/Length[normdata2];

Plot[pgram[w, normdata2, gap], {w, 0, 10}, PlotRange -> All]

enter image description here

It is clear that there are a couple of peaks and we can get refinements from NMaximize as below. I remark that it is not hard to find good bounds/starting points just by keeping values using the EvaluationMonitor option to Plot.

{m1, f1} = 
 NMaximize[{pgram[w, normdata2, gap], 2 <= w <= 4}, w]
{m2, f2} = NMaximize[{pgram[w, normdata2, gap], 4 <= w <= 6}, w]

(* Out[1434]= {2453.603870204632, {w -> 2.603084800951954}}

Out[1435]= {301.3987842825805, {w -> 5.30742917643148}} *)

Amplitudes can be estimated using a non-squared variant of the above.

amp[w_?NumberQ, vals_, gap_] := 
 2*Abs[Total[vals*Exp[I*w*gap*Range[0, Length[vals] - 1]]]]/Length[vals]

Offhand I don't know why the factor of 2 is needed in the normalization. And maybe this is not a good method anyway. But ehre is what it gives.

NMaximize[{amp[w, normdata2, gap], 2 <= w <= 4}, w]
NMaximize[{amp[w, normdata2, gap], 4 <= w <= 6}, w]

(* Out[1443]= {3.500386188934522, {w -> 2.60308480095994}}

Out[1444]= {1.226830277890204, {w -> 5.307429138313871}} *)

We even get the frequencies again.

POSTED BY: Daniel Lichtblau

I would like to make my final remark to this discussion:

But the "process" yielding the data contains exactly two vibrations, so 60-some seems a bit unsatisfactory.

It is important to notice that you do in fact not have just one ore two frequencies in your examples - as their construction suggests. There are discontinuities which lead to a whole variety of other frequencies. This is not some technical limitation of the Fourier method but a deep principle of nature.

And I think Ivan is pretty right when he says: "From my experience the main problem is to determine frequencies as accurate as possible "

Well, this is equivalent to the advice: "Measure as long as you possibly can!".

My general suggestion is: Choose one single good book about Fourier transform and its applications. There a several definitions of Fourier transform, and throughout a single book this definition is consistent - in contrast to sources in the internet! When working with Mathematica make carefully sure that you are using that same definition (by adjusting FourierParameters).

Regards -- Henrik

POSTED BY: Henrik Schachner

The code snippet below is reasonably good at figuring out how many peaks, and where they are centered in the fft.

ff2 = 0.3 + 3.5 Cos[2.6 t + .8] + 1.2 Cos[5.3 t - .4];
obstime = 80;
ffdat2 = Table[ff2, {t, 0, obstime, .1}];
fft2 = Fourier[ffdat2];
absfft2 = Abs[fft2] - Mean[Abs[fft2]];
triples = Partition[absfft2] [[1 ;; Ceiling[Length[absfft2 ]/2]]], 3, 1];
peakposns = Flatten[Position[triples, {a_, b_, c_} /; b > a && b > c]]

(* Out[1001]= {33, 68} *)

From here we want to figure out what are the actual frequencies. I take a window around each peak and weight the values in that window in order to get a finer estimate of the "true" centers. Also we subtract 1 because the DC component, at freq 0, is in position one rather than zero.

delta = 2;
centers = peakposns;
ranges = Range[centers - delta, centers + delta];
periods = Map[absfft2[[#]].(#-1)/Total[absfft2[[#]]] &, ranges]

(* Out[1005]= {32.91443987298773, 67.65933096605009} *)

We do alright but not great with this weighting.

freqs = 2*Pi*periods/obstime
Abs[{2.6, 5.3} - freqs]

(* Out[1006]= {2.58509406255003, 5.313951427743585}

Out[1007]= {0.01490593744996982, 0.01395142774358504} *)

Using a "maximal" delta of 31 gives an extremely small error for the primary frequency, around 2*10^(-5). A value of 2 gives a fairly small error for the second frequency, around 10^(-3). But I dfo not know how to assess in a methodical manner what ranges to use, and of course the above observations are based on knowing already the correct answers.

My guess is there are windowing methods better suited for this frequency recovery. In any case, as has been noted already, this is the key step. Amplitude and phase can now be recovered by linear least squares methods.

POSTED BY: Daniel Lichtblau

Hi Hans,

use the email given in the package, send $Version and your data as well.

I.M.

POSTED BY: Ivan Morozov

Hello Henrik and Ivan.

@Henrik: I fed my 2nd function in your method, and it works. But: using your really nice Manipulate-object it turns out that I need 60-some or so cosine-functions to obtain a good reproduction of the original data (assumed to be results of (error-free) measurements). But the "process" yielding the data contains exactly two vibrations, so 60-some seems a bit unsatisfactory.

Look, in my opinion the process is as follows: make your measurements or take the experimental data, apply "Fourier" and plot the spectrum. Look at the spectrum and get the number of peaks. This should be the number of periodic processes involved. In my example there are two peaks suggesting that there are two vibrations or two periodicities (as it is the case in fact). How to find them?

After all our discussion I got the impression that there is no way but to estimate frequencies (which seems to work quite well), amplitudes and phases (that works by far worse) and do a fit, and I assume yes, following Ivan, at best a linear fit with sines and cosines instead of a nonlinear fit with phases. (but I am not quite sure about this).

And I think Ivan is pretty right when he says: "From my experience the main problem is to determine frequencies as accurate as possible "

@Ivan The revised FREQUENCIES still produces lots of error messages. Of course I will send you all of them. Could you give me your email-address? And if you want to I send you the original problem as well.

Kind regards Hans

POSTED BY: Hans Dolhaine

Hi, Hans,

Yes, it may be problematic to run that package on v7, could you provide more details on the error messages? Meanwhile, are you sure you need Fourier analysis for your data? What is the process behind your data?

Perhaps, you might be interested in exploring other methods for data analysis like:

  • Wavelet analysis (ref1, ref2), but wavelet staff is available in Mathematica for v8+
  • Independent component analysis and it's variations (ref1), there is an implementation of FastICA for Mathematica somewhere on library.wolfram.com and some functions are available in Mathematica, e.g. PrincipalComponents[] and KarhunenLoeveDecomposition[] for v8+.
  • Hilbert-Huang transform (ref1), I remember seeing an implementation of the sifting process used in this method somewhere on mathematica.stackexchange.com as well as discrete Hilbert transform.
  • Dynamic Mode Decomposition (ref1), I'm not aware of any implementation in Mathematica and will be grateful if someone can provide one

This list is by no means complete, but it can be a good starting point.

P.S. I've modified the code to run (hopefully) on v7, see attached.

I.M.

Attachments:
POSTED BY: Ivan Morozov

Hello Henrik, hello Ivan,

thanks a lot to both of you.

@ Henrik: I will try it with my second example with just plugging the 2nd function in. Regarding the zeros I think you are perfectly right, that seems not to be a good idea.

@ Ivan: I managed to load FREQUENCY, but for the time being I am drowning in error messages. Maybe because I have version 7, but I will have to have a look to them.

Kind regards Hans

POSTED BY: Hans Dolhaine

Hi, Hans,

From my experience the main problem is to determine frequencies as accurate as possible and to do this usually window functions and peak interpolation is required. Take at a look at the following papers of Hunter and Laskar, where algorithms for accurate frequency identifications are described (as well as for parameters, i.e. amplitudes and phases). I personaly prefer to use Laskar's method for frequencies and the NonlinearModelFit[] for parameters, but LinearModelFit[] or LeastSquares[] can be used as well, if you search for $a_1 \cos(f_1 t) + b_1 \sin(f_1 t) + \dots$ instead of $a_1 \cos(f_1 t + \phi_1) + \dots$.

See also my attached package for signal reconstruction. It's probably not very friendly to read, but you can give it a try. You'll also need to replace CompilationTarget->"C" to CompilationTarget->"WVM" in FREQUENCY.m and remove compiler initialization.

I.M.

Attachments:
POSTED BY: Ivan Morozov

Hallo Hans,

Would you think adding zeros at the beginning and end of the data (to force it being the same at the ends of the interval) would help?

I do not think so; by adding zeros you create jumps between your data and the zeros. Maybe you add data which act as interpolation, but this is just a first (and stupid?) guess. An acquisition of a finite amount of data (on -> off) is equivalent to the application of a HeavisidePi "cut out" function. One could use a Gaussian instead. But in fact: This topic is endless! An excellent source is e.g.:

The Fourier transform and its applications /? Ronald N. Bracewell. 3rd ed. Boston : McGraw Hill, c2000.

And what about my (2nd) problem with two functions (see Fourier_MC.nb )

Well, it works exactly the same way! Just plug in the other function...

Regards -- Henrik

POSTED BY: Henrik Schachner

Hello Henrik,

that looks interesting, but (with Dan's Method) I found a reasonable frequency and amplitude in the simple case as well. The phase angle was the main problem.

Would you think adding zeros at the beginning and end of the data (to force it being the same at the ends of the interval) would help?

And what about my (2nd) problem with two functions (see Fourier_MC.nb )

Regards Hans

POSTED BY: Hans Dolhaine

Hallo Hans,

I have tried to do the same using your example.

Regards -- Henrik

Attachments:
POSTED BY: Henrik Schachner

Thanks a lot. I will have to unhurriedly have a close look at it. But what about the example in my notebook?

Would you mind to give me your email? I could send you the original problem-

Regards Hans

POSTED BY: Hans Dolhaine

Hallo Hans,

sorry for my late response! Fourier transformation is a too important topic, so your question should not be left unanswered!

I have tried to construct a hopefully meaningful example of Fourier decomposition, maybe it is helpful. Generally when dealing with Fourier it is always important to know exactly its specific meaning, and this is defined by the option FourierParameters. Have a look under "Details and Options" in the documentation.

Regards -- Henrik

Attachments:
POSTED BY: Henrik Schachner

sorry, attached the wrong notebook. take this one

Attachments:
POSTED BY: Hans Dolhaine
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract