Re: Finding simplest Fourier series between two given
- To: mathgroup at smc.vnet.net
- Subject: [mg105674] Re: [mg105662] Finding simplest Fourier series between two given
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Tue, 15 Dec 2009 07:25:34 -0500 (EST)
- References: <200912140506.AAA16985@smc.vnet.net>
Kelly Jones wrote:
> Given two Fourier series f1[x] and f2[x], where f1[x]<=f2[x] for all
> x, I want Mathematica to find the "simplest" Fourier series f3[x] that
> lies between them. More specifically:
>
> I. f1[x] <= f3[x] <= f2[x] for all x
>
> II. f3[x] has the fewest non-zero coefficients of all f3 meeting I.
>
> III. If multiple functions meet I and II, choose the one whose
> highest term is smallest (ie, the "least wiggly" one).
>
> If multiple Fourier series satisfy I, II, and III, I'll settle for any
> of them.
>
> Motivation:
>
> % I'm Fourier-fitting continuous cyclic data that's measured to the
> nearest integer. IE, a datum of 56 means the value I'm measuring is
> between 55.5 and 56.5.
>
> % Using Fourier series, I can approximate the measured data to
> arbitrary precision, but this feels silly when the terms are of order
> 0.1, 5 times smaller than the measurement precision.
>
> % I believe the data satisfies a fairly simple Fourier relation
> that's being obscured by rounding/measurement precision.
>
> "Extra credit": The data I'm measuring is cyclic, but I'm not
> necessarily feeding Mathematica an integral number or cycles. The list
> I give Mathematica may have 57.2 cycles instead of 57 or 58. The ideal
> solution would compensate for this, though I'd be happy w/ a solution
> that just works for an integral number of cycles.
I can give some code that might do what you have in mind. I cribbed part
of it from something I did several years ago, and I'm afraid I no longer
quite understand the older version, but still...
The idea is to do some thresholding to assess the most likely
approximate number of periods, remove the excess from the remaining
partial cycle, then work with a thresholded Fourier transform of that.
Here is an example. I use signal data in the range 0-255, start with
something periodic (57 cycles) but having an extra partial cycle. I then
corrupt the data muildly and by uniform random noise of mean zero. This
will not accurately reflect your situation, but may come close enough.
In[284]:= bnds = {0, 2^8 - 1};
data = RandomInteger[bnds, {20}];
new = Flatten[Table[data, {58}]];
seq = Drop[new, -14];
seq = Clip[seq + RandomInteger[{-3, 3}, Length[seq]], bnds];
Length[seq]
Out[289]= 1146
We'll look at the largest in magnitude value from the Fourier transform.
In[290]:= ft = Abs[Fourier[seq]];
Max[ft]
Out[291]= 4530.11
Now we threshold the FT so that all values less than 100 become 0, and
all values greater become 1. The point will be to assess length of runs.
In[341]:= thresh = Round[Clip[Abs[ft], {100, 100}, {0, 1}]];
lens = Map[Length, Split[thresh]]
Out[342]= {1, 54, 6, 52, 4, 55, 1, 56, 2, 55, 2, 55, 2, 56, 1, 56, 2, \
55, 2, 56, 1, 56, 2, 55, 2, 56, 1, 56, 2, 55, 2, 55, 2, 56, 1, 55, 4, \
52, 6, 54}
We pair off these runs and find the floor of the average total length.
In[346]:= newsplit = Partition[lens, 2];
cycles = Floor[Mean[Map[Total, newsplit]]]
Out[347]= 57
Now we've decided there are 57 cycles, so we remove the remaining
partical cycle from our input sequence, take the FT, and threshold it by
magnitudes to obtain mostly zeros. I set the magnitude threshold at 10
below, but 100 also seems to work well in this case. I'm not sure
offhand what is a good way to set it; maybe the cube root of the max or
some such.
In[348]:= seqnew = Take[seq, Floor[Length[seq]/cycles]*cycles];
ftnew = Map[If[Abs[#] >= 10, #, 0] &, Fourier[seqnew]];
Max[Abs[ftnew]]
Out[350]= 4532.06
As a check, we see how well we recover our sequence from this
thresholded FT.
seqapprox = Round[InverseFourier[ftnew]];
In[366]:= Max[Abs[seqapprox - seqnew]]
N[Total[Abs[seqapprox - seqnew]]/Length[seqnew]]
Length[DeleteCases[ftnew, 0]]
Out[366]= 4
Out[367]= 1.7
Out[368]= 20
The average deviation is 1.7, the max is 4, and we have 20 nonzero
frequencies. In this perhaps crude test example, we did quite well.
There is at least one other method for approximating the cycle length.
It uses lattice reduction under the hood, to find an approximate common
multiple to the binned signal values. This can work when the possible
signal range is not too large and the errors are not too large. It's
method (iii) in the psot at the URL below.
http://forums.wolfram.com/mathgroup/archive/2003/Mar/msg00184.html
This used version 4.2 or 5.0 of Mathematica, and the actual
implementation today would be a bit different.
I will mention that this does not really do what you request. It simply
tries to find a sensible period length, and delivers Fourier components
from that. You could do further thresholding to try to make more
components zero, perhaps adjust the remaining ones, and see how well you
can still approximate the signal.
Daniel Lichtblau
Wolfram Research
- References:
- Finding simplest Fourier series between two given Fourier series
- From: Kelly Jones <kelly.terry.jones@gmail.com>
- Finding simplest Fourier series between two given Fourier series