Message Boards Message Boards

1
|
4244 Views
|
1 Reply
|
5 Total Likes
View groups...
Share
Share this post:

Performing discrete Fourier transform on a TimeSeries data

Posted 2 years ago

Hi everyone, I am trying to evaluate commercial viability of wind turbines in varying weather conditions. Using the following code, I was able to get daily mean wind speeds of a coastal city in Australia for the past 10 years.

wsData = WeatherData[{"Portland", "Australia"}, 
   "MeanWindSpeed", {{2010, 1, 1}, {2020, 1, 1}, "Day"}];
DateListPlot[wsData, AspectRatio -> 1/5, Joined -> True, 
 ImageSize -> 1000, PlotStyle -> Directive[Opacity[.8], Thickness[0]]]

Daily mean wind speeds from 2010 to 2020

Now, what I would like to analyse is the underlying trend in this wind data. I would like to isolate the plot of the driving function and that of the noise in separate graphs as shown in the example I found in a text:

Discrete Fourier Transform

However, I am finding it difficult to use Fourier Transform functions with TimeSeries data, and any pointers/tips would be much appreciated.

POSTED BY: KC K

You can extract x and y values and do a bit of Fourier analysis on those.

times = wsData["Times"];
values = QuantityMagnitude[wsData["Values"]];

Check whether, or to what extent, the measurement intervals are uniform.

gaps = Differences[times];
Tally[gaps]

(* {{86400, 3575}, {172800, 27}, {432000, 2}, {259200, 3}, {345600, 1}} *)

So there were some days when the apparatus was down or measurements were missed for other reasons. Also note that the main gap, 86400, is the seconds in one day.

Now check that the missing count is what we expect.

In[45]:= Length[times]

(* Out[45]= 3609 *)

In[44]:= missingtotal = 
 Total[Cases[
   Tally[gaps], {aa_ /; aa =!= 86400, bb_} :> bb*(aa/86400 - 1)]]

(* Out[44]= 44 *)

So that will account for the full 10 years, assuming three were leap years. A clean analysis would insert values into the gaps (using the mean, perhaps) and proceed from there. We will skip this step below and do a slightly cruder analysis.

Now find the period, which, not surprisingly, is one year. We extract the position of the largest frequency from Fourier (after making the mean zero), subtract 1 to account for position one being the DC component, and divide the adjusted length by that.

len = Length[values];
truelen = len + missingtotal
ft = Abs[Fourier[values - Mean[values]]];
max = Max[ft]
mainfreq = FirstPosition[ft, max][[1]] - 1
period = N[truelen/mainfreq]

(* 3653

48.3809

10

365.3 *)

We can (again, to crude approximation) clean this by removing small frequencies. By trial and error I arrived at the factor of 0.55 below. This retains two frequencies (plus their negative mirrors, at the other end of the list).

ftModified = Map[If[Abs[#] > .55*max, #, 0.] &, ft];
Position[ftModified, aa_ /; aa > 0]

(* {{8}, {11}, {3600}, {3603}} *)

Invert this clipped frequency list and create an interpolation.

inv = Re[InverseFourier[ftModified]];
interp = Interpolation[inv];

We can plot this.

Plot[interp[t], {t, 1, len}]

enter image description here

Using the secondary frequency indicates another possible period of around 522 days, by the way.

POSTED BY: Daniel Lichtblau
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