Conference rusure::math

Check out chapter 12 of "Numerical Recipes" in your local DEC library. 

Lynn Yarbrough 

        A good text to learn applications for Fourier series would
    be an Electrical Engineering textbook on Communication Systems.
    Fourier series play an important role understanding such things
    as bandwidth, white noise, and fequency content.  The course
    I took used "Communication Systems" by Bruce Carlson.  It wasn't
    a great book but it was readable.
        I'm fairly sure that only signals which are periodic and
    have a continous first derivative have a finite number of cosine
    and sine components.  For example a square wave has an infinite
    number of harmonics, however, like most series the first few
    terms give a good approximation.
        Another problem you will have is that computing the Fourier
    series of an arbitrary signal is computationally very expensive.
    For a know signal transform tables are used.
       
    Martin H.

	Suppose you have a set of measurements (Yi,Xi), and want to
	fit these to an equation of the form:

		 N
	    Y = Sum Cj * Tj
		j=1

	Where the Tj are functions on the measured X values, and the Cj
	are to be determined by a least-squares approximation.

	Some example equations are:

	    y = c1 + c2*x + c3*x**2 + c4*x**3 + c5*x**5
	    y = c1*r*s + c2*s*t + c3*t*r + c4
	    y = c0 + c1*cos(x) + c2*cos(2*x) + c3*sin(x) + c4*sin(2*x)
	    y = c1*r*s + c2*s*t + c3*t*r + c4

	In the first case, the Tj are simply powers of the X values.
	In the third case, the Tj are functions of several (measured) values:

	    T1 = r*s, T2 = s*t, T3 = t*r, and T4 = 1.

	The problem of computing the least-squares approximations for the Cj
	means that the following expression must be minimized:

	    smp		 N
	    Sum (Y[i] - Sum Cj * Tj[i])^2
	     i		j=1

	Taking the partial derivatives w.r.t. each Cj, and setting each
	to zero, we reduce the problem to solving the following system
	of N equations (k = 1,...,N) in N unknowns (the Cjs).

	     N	       smp		   smp
	    Sum ( Cj * Sum Tk[i]*Tj[i] ) = Sum Tk[i]*Y[i],
	    j=1         i                   i

	      smp
	where Sum indicates a sum over all the sample data.
	       i

    I'll second the recommendation for a peek at "Recipes" for a few
    simple Fourier analysis routines.

    Good information on digital signal procesing is also available in
    Bracewell's book on Fourier integrals, though not with an emphasis
    on computational (programming) aspects.  I can recommend many other
    books on signal processing if necesary.  The subject is very easy
    and concrete, so should be a snap for anyone with 'mathematical
    maturity'.

    Re .-1 - though that is a neat way of looking at the problem from
    first principles, it neglects the crucial benefits of (complete)
    orthogonal basis functions;  these are what make the Fourier
    technique and its relatives so valuable.

    Your problem sounds like classical time series analysis - what is
    it more specifically?  I may be able to make a more accurate
    recommendation on a line of attack.

    - Jim

	If you have a continous or regular sampling of points,
	and the Tjs are orthogonal (over the range of sample points)
	then the term:

	      	       smp
	               Sum Tk[i]*Tj[i]
	                i

	can be greatly simplified.  If the Tjs are orthonormal, then
	it's even better.

    Try the code in 760.2, it is easy to use:
    
    Take your real valued data array and put it in the real "FR" array.
    
    Zero out the imaginary "FI" array.
    
    Call the routine, with the proper length "N" parameter, i.e. 8 for
    256 points, for example.
    
    You get back a real and imaginary result. Find the square root of
    the sum of the squares of the real and imaginary array elements.
    
    i.e.  ans(i)=sqrt(fr(i)**2+fi(i)**2)
    
    You will find that the result array is symmetric from the lower
    half to the upper half of the  array, so you only have to
    use the lower half. 
                      
    The values you get will represent the amplitude at each succesive
    frequency component. This data masks the relative phase difference
    between each frequency component, but for many analysis efforts
    the phase of each frequency (or "time starting point") is not really
    that important. For communications problems the phase is often
    critical.
     
    You will also find that that result often has a wide dynamic range,
    with values ranging over five or six orders of magnitude. It is
    common to take the base-10 log of the data, for display purposes,
    then the results will range more narrowly. This does not distort
    the intuitive results you are looking for.          
    

Thanks for the contributions so far.  Here's what I've done.  Suppose we
have the following 16 data points, representing weekly sunspot counts for
the past 16 weeks.

			  **
	  40	  *      *  *
		 * *    *
	  38	*   *  *       *
		     **	     **
	  36		      
		----+----+----+---
		    5   10   15
		   Week number

I'd like to know what future week has the best chance of catching the most 
sunspots.  It might be 3, 4, or 5 weeks from today I would guess, from eyeing
the data.  But I can't fly to Chile (OK, Arizona) and wait around 3 weeks
and pick the best set of data -- solar telescope time has to be scheduled.

So taking the data and feeding it into the FFT routine in "Numerical 
Recipes",  I get out a complex array.  As suggested in .6, I compute the value
of the "ans" array, and it does look rather symmetric.  Here's the output:


 Index  Input   Trans-R     Trans-I     "ans(i)"  

     1  38.00   618.00000     0.00000   618.00000
     2  39.00    -3.39635    -3.93755     5.19995
     3  40.00     0.70711    13.36396    13.38266
     4  39.00     1.17218     0.13438     1.17986
     5  38.00    -1.00000    -3.00000     3.16228
     6  37.00    -1.75797    -0.45141     1.81500
     7  37.00    -0.70711    -0.63604     0.95108
     8  38.00    -0.01786    -0.52334     0.52364
     9  39.00     0.00000     0.00000     0.00000
    10  40.00    -0.01786     0.52334     0.52364
    11  41.00    -0.70711     0.63604     0.95108
    12  41.00    -1.75797     0.45141     1.81500
    13  39.00    -1.00000     3.00000     3.16228
    14  37.00     1.17218    -0.13438     1.17986
    15  37.00     0.70711   -13.36396    13.38266
    16  38.00    -3.39635     3.93755     5.19995

OK, folks, where do we go from here?  How do I calculate "projected"
values for weeks 17 and beyond?  Where are these sine and cosine
coefficients I'm supposed to get out?  Am I in left field?

  John

    Your data is probably a bit sketchy to make really confident
    predections.  But essentially, if you analyze some time series and
    it shows periodic behaviour, (that is, some sharp peaks in the
    frequency domain),  you can try extrapolating those spectral
    components beyond the time interval you started with and (hopefully)
    predict the future.

    However, I believe sunspot activity has a large 1/f noise component
    which would make such predictions rather suspect.

    Do you possibly have any more extensive collections of data?  Then
    you could try the above suggestion, and see if some early parts of
    the data do succeed in predecting the tail end of the data set.

    - Jim

    Can't we just solve my problem -- learning -- and not complain about
    the data source?  For educational purposes it would seem that 16
    points is a lot easier to deal with than, say, 256 or 8192.
    
>    ... you can try extrapolating those spectral
>    components beyond the time interval you started with and (hopefully)
>    predict the future.

    Well it would appear that the data in .6, sketchy though it may be,
    shows periodic behavior.  So now what does this Fourier analysis
    buy me?  Realizing full well that the DATA is insufficient, what
    is the METHOD one uses?  What does it mean to "extrapolate those
    spectral components"?  Whence cometh coefficients Ai, Bi, Ci, Di,
    (or whatever) for the best fit of the form:
    
    	16?
    Y =	Sum  Ai * SIN(Bi * X) + Ci * COS(Di * X)
    	i=1

    Is this the wrong question?  Wouldn't one normally construct such
    an equation to predict future values from periodic prior points?
    Isn't that what Fourier analysis is all about?
    
      John

    I didn't mean to 'complain' about the data source, just mention that
    it may not be enough to get a very reliable esitmate.

    What you probably want to try is linear prediction, which is also briefly
    covered in the Recipes book.  Make the assumption that the data results
    from the impulse response of some undertermined infinite impulse response
    digital filter (also called a recursive filter).

    One example of such a response is the Fibonacci sequence.

    The idea is to assume your filter has a certain number of poles in the
    complex plane and fit those poles, and the initial conditions, to a
    supply of data.  Then let the filter run for a while longer to extrapolate
    into the future.

    The Recipes book gives routines for doing maximum entropy estimation and
    then using that in a linear prediction process, so it should be pretty
    much a plug in procedure.

    Fourier series (such as the DFT) makes the assumption that the data
    in the 'window' is periodic; that is the content of the window repeats
    periodically outside the supply of data.  If the data at the endpoints
    of the window do not smoothly join, (such as occurs when a sine wave
    has a non-integral number of cycles fitting in the window), there
    will be jagged discontinuities in the periodic time domain signal that
    the DFT is assuming.  The idea of the various windows, such as Hamming,
    is to force the periodic repetition of the data to be joined smoothly
    by gradually attenuating the data near the endpoints, and cut down on
    the harmonics that result.  The problem with the DFT is that very low
    frequencies, such as only one or two cycles, cannot be accurately
    resolved because the DFT samples the signal in both the time and
    fequency domains;  a very low time domain spectral component will 'alias'
    and fall between the samples in the frequency domain.  DFT only works
    well when many cycles of a periodic spectral component are available
    and cannot resolve fractional cycle frequencies at all.

    This is the problem you face with the sunspot data, so that the DFT
    is not really the correct tool.  I think the maximum entropy technique
    is far better for this problem.  One possible fudge to try is oversample
    your data supply by fitting a spline curve to the data points and then
    feed that data (with additional points interpolated on the spline) to
    the spectral analysis.  It might increase the accuracy of the esitmate
    (but then it may not!)

    But again - with such a small amount of data you could probably eyeball
    the graph and estimate the future as well or better than a computer
    program can!  Most of these magical procedures cannot really work
    miracles.  It's widely acknowledged that spectral esitmation is really
    an art as much as a science in the real world...

    - Jim

    I tried feeding the data above into a maximum entropy estimation
    routine and then letting the resulting recurrance run for about 10
    more sample points;  varying the number of coefficients (poles)
    of the filter over the range 4 to 10 gave fairly stable estimates
    of peak sunspot activity at 3 weeks into the future, with a minimum
    just toward the end of the 6'th week.

    It will be interesting to check the predictions against what actually
    happens...

    - Jim

    Could you let me know when to drop out of the stock market?  Thanks.

>    Could you let me know when to drop out of the stock market?  Thanks.

    Close of trading, August 25th, 1987.