# Performing discrete Fourier transform on a TimeSeries data

Posted 1 year ago
1974 Views
|
|
5 Total Likes
|
 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]]] 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:However, I am finding it difficult to use Fourier Transform functions with TimeSeries data, and any pointers/tips would be much appreciated.
 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}] Using the secondary frequency indicates another possible period of around 522 days, by the way.