1
$\begingroup$

I would like to extract the most important frequency modes from a data set which exhibits strong annual periodicity, as well as some (less) important shorter-term frequency components. The following approximates the data series fairly well:

n = 3*365;
data = If[# < 0.3, 0, #] & /@ MovingAverage[1.6 #^2 & /@ (Table[Sin[ \[Pi] x/365], {x, n}] + 
   RandomReal[{0.5, -0.5}, {n}]), 20];
ListLinePlot[data]

Generating:

enter image description here

I have tried taking the discrete Fourier transform of the data, and then looking at its absolute value to identify the most important frequency components using:

ListPlot[Abs[Fourier[data]]]

However, I am having trouble 'seeing', then extracting, the most important components, which I hope should correspond to the annual periodicity, as well as some shorter-term components, of perhaps a few months in period.

Does anyone know how I can do this? I know this probably as much to do with my lack of knowledge about Fourier series, as it is about my Mathematica inexperience!

Update: to be clear, I would like to extract the most important modes, then reconstruct the data series using only these few.

Best,

Ben

$\endgroup$
1
  • $\begingroup$ Try Periodogram. $\endgroup$
    – bbgodfrey
    Feb 24, 2015 at 13:26

1 Answer 1

3
$\begingroup$

To use Fourier, subtract out the mean (otherwise the "DC" component swamps out the period you are trying to see).

q = data - Mean[data];
ListLinePlot[Abs[Fourier[q]][[1 ;; 20]], PlotRange -> All]

enter image description here

I only plotted the first few terms so it would be easier to see the low frequency region. What you see here is a large peak in the 4th bin (the first is DC, the second is one oscillation across the data set, etc.) So the peak at 4 represents a fundamental periodicity of 3 oscillations across the data set. Since your data set is 3*365 long, this is one oscillation per 365 samples.

If you wish to go a little deeper, you can reconstruct the signal using variable numbers of Fourier parameters. Here is a Manipulate that demonstrates (make sure to define the data sequence first). The heart of this is the line rec which builds an approximation to the data signal using the magnitude and phase of the Fourier coefficients and the Cos at the appropriate frequencies.

Manipulate[Module[{}, sig = data - Mean[data]; lenSig = Length[sig];
  fftY = Fourier[sig, FourierParameters -> {1, -1}];
  mag = Abs[fftY]/lenSig; phase = Arg[fftY];
  ordMag = Ordering[mag, All, GreaterEqual];
  rec = Total[Table[mag[[ordMag[[k]]]] Cos[2 Pi (ordMag[[k]] - 1) (n - 1)/lenSig + 
        phase[[ordMag[[k]]]]], {k, 1, numTerms}, {n, 1, lenSig}]];
  GraphicsRow[{ListPlot[{Tooltip[mag, "not used"], 
      Tooltip[Thread[{ordMag[[1 ;; numTerms]], 
         mag[[ordMag[[1 ;; numTerms]]]]}], "used in reconstruction"]},
      Filling -> Axis, PlotRange -> All, PlotLabel -> "Fourier Transform of Signal", 
      PlotStyle -> {Blue, {Black, PointSize[0.015]}}], 
    ListLinePlot[{Tooltip[sig, "signal"], 
      Tooltip[rec, "reconstruction"]}, 
      PlotLabel -> "Signal (blue) and Reconstructed Signal (brown)", 
      PlotRange -> All, Filling -> {1 -> {2}}]}, ImageSize -> 800]], 
   {{numTerms, 1, "number of terms in reconstruction"}, 1, 30, 1, Appearance -> "Labeled"}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Great - thanks for that. Can I ask how I can then go about extracting the Fourier coefficients, for the most important few modes, then reconstructing the series using these? Sorry for my lack of knowledge here! $\endgroup$
    – ben18785
    Feb 24, 2015 at 14:06
  • $\begingroup$ That's great; many thanks for your help here. I really appreciate the time taken for you to make such an in-depth response; particularly one that allows me to dynamically look at the answer. Best, Ben $\endgroup$
    – ben18785
    Feb 24, 2015 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.