How to find the frequencies of a time series which contain datetime ?

How to find the frequencies of a time series which contain datetime, as in the following example ?

(Note: I did not find such an exmple in the Fast Fourier transform (fft) page)

a = {
    '17-Jun-2021 12:00:00', 25
    '18-Jun-2021 20:00:00', 21
    '19-Jun-2021 06:00:00', 31
    '20-Jun-2021 14:00:00', 35
    '21-Jun-2021 17:00:00', 31
    '22-Jun-2021 16:00:00', 41
    '23-Jun-2021 16:00:00', 22
    '24-Jun-2021 13:00:00', 22
    '25-Jun-2021 15:00:00', 45
    '26-Jun-2021 11:00:00', 20
    '27-Jun-2021 19:00:00', 25
    '28-Jun-2021 18:00:00', 25
    '29-Jun-2021 15:00:00', 29
    '30-Jun-2021 10:00:00', 31
    '01-Jul-2021 15:00:00', 20
    '02-Jul-2021 10:00:00', 28
    '03-Jul-2021 14:00:00', 31
    '04-Jul-2021 13:00:00', 24
    '05-Jul-2021 09:00:00', 24
    '06-Jul-2021 14:00:00', 25
    '07-Jul-2021 15:00:00', 40}
plot(datetime(a(:,1)),cell2mat(a(:,2)))

Answer 1

[Abderrahim. B]

Hi!

You can use pspectrum function. It returns both power and frequency arrays.

% Your data
a = {
    '17-Jun-2021 12:00:00', 25
    '18-Jun-2021 20:00:00', 21
    '19-Jun-2021 06:00:00', 31
    '20-Jun-2021 14:00:00', 35
    '21-Jun-2021 17:00:00', 31
    '22-Jun-2021 16:00:00', 41
    '23-Jun-2021 16:00:00', 22
    '24-Jun-2021 13:00:00', 22
    '25-Jun-2021 15:00:00', 45
    '26-Jun-2021 11:00:00', 20
    '27-Jun-2021 19:00:00', 25
    '28-Jun-2021 18:00:00', 25
    '29-Jun-2021 15:00:00', 29
    '30-Jun-2021 10:00:00', 31
    '01-Jul-2021 15:00:00', 20
    '02-Jul-2021 10:00:00', 28
    '03-Jul-2021 14:00:00', 31
    '04-Jul-2021 13:00:00', 24
    '05-Jul-2021 09:00:00', 24
    '06-Jul-2021 14:00:00', 25
    '07-Jul-2021 15:00:00', 40}

a = 21×2 cell array
    {'17-Jun-2021 12:00:00'}    {[25]}
    {'18-Jun-2021 20:00:00'}    {[21]}
    {'19-Jun-2021 06:00:00'}    {[31]}
    {'20-Jun-2021 14:00:00'}    {[35]}
    {'21-Jun-2021 17:00:00'}    {[31]}
    {'22-Jun-2021 16:00:00'}    {[41]}
    {'23-Jun-2021 16:00:00'}    {[22]}
    {'24-Jun-2021 13:00:00'}    {[22]}
    {'25-Jun-2021 15:00:00'}    {[45]}
    {'26-Jun-2021 11:00:00'}    {[20]}
    {'27-Jun-2021 19:00:00'}    {[25]}
    {'28-Jun-2021 18:00:00'}    {[25]}
    {'29-Jun-2021 15:00:00'}    {[29]}
    {'30-Jun-2021 10:00:00'}    {[31]}
    {'01-Jul-2021 15:00:00'}    {[20]}
    {'02-Jul-2021 10:00:00'}    {[28]}

% Your plot
plot(datetime(a(:,1)),cell2mat(a(:,2)))

% convert to table, then to timetable
aTbl = cell2table(a) ;
aTbl.Properties.VariableNames = ["Time" , "timeSerie"] ;
aTbl.Time = datetime(aTbl.Time, "InputFormat","dd-MMM-uuuu HH:mm:ss") ;
aTimeTbl = table2timetable(aTbl) ;

% Some preprocessing
aTimeTBL_RETIME = retime(aTimeTbl, "hourly", "previous") ;

% spectrum : returns power and frequency vectors
[pow, freqs] = pspectrum(aTimeTBL_RETIME) 

pow = 4096×1

1.0e+03 * 1.4474 1.4449 1.4373 1.4249 1.4077 1.3858 1.3595 1.3291 1.2947 1.2568

freqs = 4096×1

1.0e-03 * 0 0.0000 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0003 0.0003

% Plot power in log scale vs frequency
figure
semilogy(freqs, pow)

Hope this helps

[Sim]

@Issa thanks!

Would you be able to tell me what is the dominant frequency in the last plot ?

(that plot is not clear to me..)

semilogy(freqs, pow)

[dpb]

The "dominant" frequency is the DC component.

As noted, you really can't find one and the overall shape is basically just the typical 1/F rolloff.

The attached shows what was taken the PSD of -- it's nothing but a bunch of step functions; the blue line is the interpolated/retimed signal, the red asterisks are the only data points there actually are -- and trying to interpolate something over a period of two weeks sampled sporadically at maybe once a day is a fool's errand.

NB: that removing the mean does eliminate the DC so one can see what there is there, but again, I think it is meaningless.

"Just because you CAN do something doesn't mean you should."

ADDENDUM/ERRATUM:

I had inadvertently scaled the frequency axis before by 3600 sec/hr which screwed up the frequencies -- with such short sample length but such a slow sample rate, Fmax = 1/(2dt) --> 1/(2*24.15) --> 0.0207 hr-1, not the 0.5 hr-1 that was shown before.

[Sim]

Thanks a lot @dpb!

About "The "dominant" frequency is the DC component."

Let me write a Fourier series here:

As far as I know or I remember, the so-called “Direct Current” (DC) component is used to indicate the

term in the Fourier series (probably in electrical engineering?), i.e. the “zero frequency component”.

In the frequency-domain representation, i.e. your second plot, the DC component would be then:

%  Magnitude ~ 1 (or maybe ~ 2), F = 0 cycles/hour

If I still remember correctly, the so called "dominant frequency" is the frequency with the largest amplitude (magnitude), i.e. that one you have indicated in your second plot:

% Magnitude ~ 19.21, F = 0.017 cycles/hour

OK, so if the dominant frequency is F = 0.017 cycles/hour, then the period of one cycle would be

% T = 1 / F  = 58.8 hours ~  2.5 days

About the "As noted, you really can't find one and the overall shape is basically just the typical 1/F rolloff."

I did not understand what you meant...

About the "The attached shows what was taken the PSD of ..."

Could you please add the code you have used to create your graphs as well ?

About the "NB: that removing the mean does eliminate the DC so one can see what there is there, but again, I think it is meaningless."

I did not understand what you meant...

[Abderrahim. B]

Hi! semilogy( X , Y ) plots x- and y-coordinates using a linear scale on the x-axis and a base-10 logarithmic scale on the y-axis. I used it because power usually plotted in logarithmic scale. You can just use normal plot if you semilogy not relevant to what you want.

[Sim]

Thanks for the addendum @dpb!

Just for my curiosity and to understand what you did, would you be so kind to write here the code for this plot please ?

Did you use fft or pspectrum ?

Sorry, but I am getting confused, and I did not understand your reasoning....

From your plot I read that the largest amplitude of 19.21 occures at the dominant frequency of

% F = 0.00069

With such a dominant frequency, the period would be

% T = 1 / 0.00069 ~ 1449 seconds = 0.4 hours ~ 24 minutes

But then you wrote also, in your previous comment that:

% Fmax = 1/(2dt) --> 1/(2*24.15) --> 0.0207 hr-1

which would give a period of

% T = 1 / 0.0207 ~ 48 hours

Then, is Fmax the dominant frequency ? or something different ?

Probably with a code attached here I could understand better.... :-)

[Sim]

Thanks @Issa for clarifying :-)

[dpb]

"About "The "dominant" frequency is the DC component."

All I was doing was pointing out that if you don't detrend/demean the input signal the DC component is by far the largest compoenent in the resulting PSD -- and it clouds looking at the plot when what little structure there is is in the very nearby neighborhood. No more, no less...and you'll note I put "dominant" in quotes to emphasize the point that it was "just the DC" and not really a peak.

But, yes, one can do the same thing as subtracting the mean from the time series by simply zero'ing out the DC bin in the resulting PSD/FFT; it has no effect on the computation -- just like each frequency bin is not dependent upon the magnitude of the adjacent ones.

You have to admit that removing it (DC) from the plot makes the rest much more easily visible -- in particular, what it does is then connect the origin point to the first bin with a line back to the origin that visually emphasizes that first lobe.

I totally misunderstood in part that these weren't the only data you might have so I did get focused on the fact that only 20 points is simply an insufficient number of samples from which to even begin to try to do spectral analysis -- but the signal still has to be terribly undersampled for any variable that I can think one would be measuring -- what is this data measuring that would only change once a day??? -- and it's pretty clear just looking at what data you do have that it does change quite a bit from day to day; does it really jump from 22 on day 7 to 45 on day 8 instantaneously (or nearly so) on day 8? One would not think that to be very likely behavior; and one would also think quite probably there are other trends going on as well.

Hence, the conclusion that it is probably undersampled drastically follows -- that may or may not make any difference depending upon what the data are and what conclusions are being attempted to be drawn from it -- that we know nothing about.

As far as the Nyquist frequency, it is solely dependent upon the actual sampling frequency -- it is defined as 1/2 the sample frequency and has nothing to do with the real signal at all as far as defining what it is for a given dataset. Where it does come into play and what it means is you canNOT estimate any components of any signal component that is also in the raw data that you're sampling that is any higher in its frequency content than what the Nyquist frequency is -- and those components, if they exist in the raw data, are aliased back into the computed FFT and cannot be distinguished from nor separated from the real input in the same (aliased) frequency bin.

This means, simply put, that if there were any fluctuations between the collected daily points you have of the variable, then those components HAVE been included in that data but are not able to be separated out into their rightful higher-frequency locations.

So, my concern still exists that even if you have 10X (or 100X) the number of points that without a lot more information about what you're data are to represent that you need to be very cautious about drawing any real conclusions from these data.

Frequency components above the Nyquist frequency can ONLY be removed from the data BEFORE it is collected and by anti-aliasing filters in the analog domain -- once the data are sampled, IT'S TOO LATE AND YOU CAN'T FIX IT.

<Wikipedia Aliasing> is not a too bad intro with some examples in audio/visual as well as graphics -- and it'll be covered in any good text of digital signal processing; I'm such an old fogey all my favorites are 40-yr out of print so I'll not even try something modern.

[Sim]

Thanks a lot for your articulated answer @dpb, very kind!

I was just exploring one potential analysis in this direction (which looks like a typical tool for electrical engineering), but I understood that these data are simply not enough to perform a meaningful analysis.. I will think about something else then - Very grateful for this positive critique! :-)

[dpb]

Oh, why didn't you say so!!! :) Wouldn't have gotten so uptight on trying to prevent a perhaps egregious blunder... :)

OK, for pedgogical purposes, I'd suggest to look at the sample datasets with MATLAB used in the documentation as starters...there's a sunspot dataset that looks at longer-time calendar data -- although it's not in datetime format (its inclusion having preceded its introduction by 30 years), you could convert the yearly data to datetimes and to a timetable if you wanted practice using the newfangled datetime and other features as well...although truthfully for spectral and time-series analysis, the reversion back to plain doubles is often more convenient because all the signal processing routines haven't been converted to know anything about datetimes or durations. timetable has, but not necessarily the core Signal Processing routines.

As far as figuring out the frequency of the low-sampling frequency signals, to see have it right you can always convert the times to seconds and deal in Hz and then scale the resulting frequencies back to the longer time frames afterwards. It is easy to get confused with days and years in those locations instead, granted.

But, remember and commit to memory the base fundamental relationships in data sampling/signal analysis -- if you remember these, you can always get it right --

Fs = 1/dt                - Sampling frequency is 1/delta-t
Fmax = 1/(2*dt) (= Fs/2) - Fmax (Nyquist) is one-half sampling rate
T = N*dt (= 1/df)        - Sample time span is number samples * dt
df = Fmax/(N/2)          - frequency resolution is Fmax over half the number samples

Applying those will get you out of the morass; if one of them doesn't hold you've made a mistake somewhere!

[Sim]

thanks a lot @dpb for your explanation, very appreciated! :-)

Answer 2

[dpb]

I think it's absurd to even think about peforming spectral analysis on such a intermittently-sampled signal, but the only thing even close to a dominant frequency would be at roughly 0.015/hr and could only be distinguished at all after detrending the trace to remove DC component.

To believe this would be of any real meaning, however, would be nuts...

[Sim]

Thanks for your comment @dpb, very kind!

About "I think it's absrud to even think about peforming spectral analysis on such a intermittently-sampled signal"

Do you mean that the signal is undersampled, i.e. there are not enough data to compose that signal, which would not satisfy the Nyquist theorem's criterion ?

[dpb]

There simply isn't enough data at all -- and what sampling there is is hit-n-miss as far as timing -- and we have no idea what the sampled data are, but guessing something for a once-a-day value to interpolate in the values is nothing more than that -- just a guess unless there's somethng known about what it is that is being recorded that makes it even feasible.

>> tD.Date.Format='dd-MMM-uuuu HH'
tD =
  21×4 table
         Date         V          T           dT  
    ______________    __    ___________    ______
    17-Jun-2021 12    25         0 days     32 hr
    18-Jun-2021 20    21    1.3333 days     10 hr
    19-Jun-2021 06    31      1.75 days     32 hr
    20-Jun-2021 14    35    3.0833 days     27 hr
    21-Jun-2021 17    31    4.2083 days     23 hr
    22-Jun-2021 16    41    5.1667 days     24 hr
    23-Jun-2021 16    22    6.1667 days     21 hr
    24-Jun-2021 13    22    7.0417 days     26 hr
    25-Jun-2021 15    45     8.125 days     20 hr
    26-Jun-2021 11    20    8.9583 days     32 hr
    27-Jun-2021 19    25    10.292 days     23 hr
    28-Jun-2021 18    25     11.25 days     21 hr
    29-Jun-2021 15    29    12.125 days     19 hr
    30-Jun-2021 10    31    12.917 days     29 hr
    01-Jul-2021 15    20    14.125 days     19 hr
    02-Jul-2021 10    28    14.917 days     28 hr
    03-Jul-2021 14    31    16.083 days     23 hr
    04-Jul-2021 13    24    17.042 days     20 hr
    05-Jul-2021 09    24    17.875 days     29 hr
    06-Jul-2021 14    25    19.083 days     25 hr
    07-Jul-2021 15    40    20.125 days    NaN hr
>> 

is your data with a column of the time from first sample between points, T and then the dT of those -- it goes from a low of 10 hrs to as high as 32 hrs between samples.

The average sampling rate is 24.15 hr.

The more nearly what one can get from such a paucity of data looks more like --

The blue is the actual data itself -- very little in the way of structure could be claimed, the yellow-gold(?) is the result of interpolating in the frequency domain by augmenting the 21 data points with enough zeros to compute a 1024-pt FFT. Introduces a lot of spurious peaks; looks nice and smooth but you couldn't believe any of it is really real. The maroon is from the NUFFT() using the actual sampled times; it's only maringally different in overall result from the discrete although one might claim there's a peak around 0.006 hr-1 -- but, again, you'd be kidding yourself if you beliveed it was real.

There's simply not enough data to even think of this...

[Sim]

Thanks a lot @dpb, clear.. :-) Actually, I would have more data but I did not post them here in this example.. :-)

P.S.: Still for clarity, putting the code or the function names to get the third and fourth columns, i.e. T and dT, in your table, would be nice.

For examples, the fourth column, i.e. dT could be derived in this way:

>> hours(diff(datetime(a(:,1))))

ans =

    32
    10
    32
    27
    23
    24
    21
    26
    20
    32
    23
    21
    19
    29
    19
    28
    23
    20
    29
    25

[dpb]

The code to generate dT actually was more like

tD.T=tD.Date-tD.Date(1);    % the time vector
tD.T.Format='h';            % set the format
tD.dT=[diff(tD.T);nan];     % must have height(tD) elements in table

Salt to suit...

[Sim]

thanks :)