convert time domain data to frequency domain

Dear all,

I want to convert this data (time-domain) to frequency domain using plomb and fft. The structure of data I have like this:

Col 1: Year, Col 2: Month, Col 3: Date, Col 4: Hours, Col 5: Data

The interval between data is 1 hour, even though there are missing data. How can I convert it to frequency domain using plomb or fft?

Thank you!

[Adib Muhammad]

the frequency must be in cycle per day

Answer 1

[Star Strider]

There are 6839 missing (non-sequential) rows in this file, so to use fft you will first have to interpolate them so they are sequentiial and have a constant (hourly) sampling interval. Yoiu can do that using table2timetable and then retime. The plomb function does not need consatant sampling intervals, and can use datetime arrays for its time variable, so that is straightforward. The fft calculation requires a bit of pre-processing.

This is a difficult problem reequiring that the data in the fiirst four columns be converted to a datetime array, so I went ahead and did all of it.

If you want to calculate the Fourier transform on the original (non-resampled) data in ‘T1’ or ‘TT1’ (not ‘TT1r’), use the nufft function. II leave that to you.

T1 = readtable("Time_domain_data.xlsx");
T1.Properties.VariableNames = {'Year','Month','Day','Hour','Data'}

T1 = 136755x5 table
    Year    Month    Day    Hour     Data 
    ____    _____    ___    ____    ______

    2004      1       1       1     2.6066
    2004      1       1       2      2.602
    2004      1       1       3     2.5964
    2004      1       1       4     2.6042
    2004      1       1       5     2.6075
    2004      1       1       7     2.6224
    2004      1       1       8     2.6254
    2004      1       1       9      2.626
    2004      1       1      10     2.6165
    2004      1       1      11     2.6131
    2004      1       1      12     2.5956
    2004      1       1      13      2.584
    2004      1       1      14     2.5718
    2004      1       1      16     2.5286
    2004      1       1      17     2.5223
    2004      1       1      18     2.5122

% yu = unique(T1{:,1})
% mu = unique(T1{:,2})
% du = unique(T1{:,3})
% hu = unique(T1{:,4})
DateTime = datetime([T1{:,1:4} zeros(size(T1,1),2)], Format="yyyy-MMM-dd HH");
T1 = addvars(T1,DateTime, 'Before',1);
T1 = removevars(T1,2:5)

T1 = 136755x2 table
       DateTime        Data 
    ______________    ______

    2004-Jan-01 01    2.6066
    2004-Jan-01 02     2.602
    2004-Jan-01 03    2.5964
    2004-Jan-01 04    2.6042
    2004-Jan-01 05    2.6075
    2004-Jan-01 07    2.6224
    2004-Jan-01 08    2.6254
    2004-Jan-01 09     2.626
    2004-Jan-01 10    2.6165
    2004-Jan-01 11    2.6131
    2004-Jan-01 12    2.5956
    2004-Jan-01 13     2.584
    2004-Jan-01 14    2.5718
    2004-Jan-01 16    2.5286
    2004-Jan-01 17    2.5223
    2004-Jan-01 18    2.5122

dt = diff(hour(T1{:,1}));
nonseq = nnz(dt>1)                                                  % Number Of Non-Sequantial Times

nonseq = 6839

figure
plot(T1.DateTime, T1.Data)
grid
xlabel('Date & Time')
ylabel('Data')
title('Table Of Original Data')

[pxx,freq] = plomb(T1.Data, T1.DateTime);
[pxxmax,idx] = max(pxx);
fprintf('\nThe PSD of the peak at %.3E 1/hour is %.3f dB/Hz and the period is %.3E hours\n\n',freq(idx), pxxmax, 1/freq(idx))

The PSD of the peak at 3.169E-08 1/hour is 521346.450 dB/Hz and the period is 3.156E+07 hours

figure
semilogx(freq, pxx)
grid
xlabel('Frequency (Hour^{-1})')
ylabel('Power / Frequency  (dB/Hz)')
title('''plomb''')

TT1 = table2timetable(T1);
TT1r = retime(TT1,'hourly','pchip');
s = TT1r.Data;
t = 0:size(TT1r,1)-1;                                               % Cumulative Hours

[FTs1,Fv] = FFT1(s,t);

[magmax,idx] = max(abs(FTs1)*2);
fprintf('\nThe magniitude of the peak at %.3E 1/hour is %.3f and the period is %.3E hours\n\n',Fv(idx), magmax, 1/Fv(idx))

The magniitude of the peak at 1.144E-04 1/hour is 0.041 and the period is 8.738E+03 hours

figure
semilogx(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hour^{-1})')
ylabel('Magnitude)')
title('''fft''')

% % % % % ONE-SIDED FOURIER TRANSFORM — 
function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

.

[Adib Muhammad]

thank you for your answer, but I'm curious, the frequency is "frequency (hour-1)". Is that possible to make the frequency in 'cycle per day'?

[Star Strider]

My pleasure!

Multiply the frequency axis vector by 24. There are more cycles per day than per hour, so the frequency of the peak increases proportionally.

T1 = readtable("Time_domain_data.xlsx");
T1.Properties.VariableNames = {'Year','Month','Day','Hour','Data'}

T1 = 136755x5 table
    Year    Month    Day    Hour     Data 
    ____    _____    ___    ____    ______

    2004      1       1       1     2.6066
    2004      1       1       2      2.602
    2004      1       1       3     2.5964
    2004      1       1       4     2.6042
    2004      1       1       5     2.6075
    2004      1       1       7     2.6224
    2004      1       1       8     2.6254
    2004      1       1       9      2.626
    2004      1       1      10     2.6165
    2004      1       1      11     2.6131
    2004      1       1      12     2.5956
    2004      1       1      13      2.584
    2004      1       1      14     2.5718
    2004      1       1      16     2.5286
    2004      1       1      17     2.5223
    2004      1       1      18     2.5122

% yu = unique(T1{:,1})
% mu = unique(T1{:,2})
% du = unique(T1{:,3})
% hu = unique(T1{:,4})
DateTime = datetime([T1{:,1:4} zeros(size(T1,1),2)], Format="yyyy-MMM-dd HH");
T1 = addvars(T1,DateTime, 'Before',1);
T1 = removevars(T1,2:5)

T1 = 136755x2 table
       DateTime        Data 
    ______________    ______

    2004-Jan-01 01    2.6066
    2004-Jan-01 02     2.602
    2004-Jan-01 03    2.5964
    2004-Jan-01 04    2.6042
    2004-Jan-01 05    2.6075
    2004-Jan-01 07    2.6224
    2004-Jan-01 08    2.6254
    2004-Jan-01 09     2.626
    2004-Jan-01 10    2.6165
    2004-Jan-01 11    2.6131
    2004-Jan-01 12    2.5956
    2004-Jan-01 13     2.584
    2004-Jan-01 14    2.5718
    2004-Jan-01 16    2.5286
    2004-Jan-01 17    2.5223
    2004-Jan-01 18    2.5122

dt = diff(hour(T1{:,1}));
nonseq = nnz(dt>1)                                                  % Number Of Non-Sequantial Times

nonseq = 6839

figure
plot(T1.DateTime, T1.Data)
grid
xlabel('Date & Time')
ylabel('Data')
title('Table Of Original Data')

[pxx,freq] = plomb(T1.Data, T1.DateTime);
[pxxmax,idx] = max(pxx);
fprintf('\nThe PSD of the peak at %.3E 1/hour is %.3f dB/Hz and the period is %.3E hours\n\n',freq(idx), pxxmax, 1/freq(idx))

The PSD of the peak at 3.169E-08 1/hour is 521346.450 dB/Hz and the period is 3.156E+07 hours

figure
semilogx(freq, pxx)
grid
xlabel('Frequency (Hour^{-1})')
ylabel('Power / Frequency  (dB/Hz)')
title('''plomb''')

figure
semilogx(freq*24, pxx)
grid
xlabel('Frequency (Day^{-1})')
ylabel('Power / Frequency  (dB/Hz)')
title('''plomb''')

TT1 = table2timetable(T1);
TT1r = retime(TT1,'hourly','pchip');
s = TT1r.Data;
t = 0:size(TT1r,1)-1;                                               % Cumulative Hours

[FTs1,Fv] = FFT1(s,t);

[magmax,idx] = max(abs(FTs1)*2);
fprintf('\nThe magniitude of the peak at %.3E 1/hour is %.3f and the period is %.3E hours\n\n',Fv(idx), magmax, 1/Fv(idx))

The magniitude of the peak at 1.144E-04 1/hour is 0.041 and the period is 8.738E+03 hours

figure
semilogx(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hour^{-1})')
ylabel('Magnitude)')
title('''fft''')

figure
semilogx(Fv*24, abs(FTs1)*2)
grid
xlabel('Frequency (Day^{-1})')
ylabel('Magnitude)')
title('''fft''')

% % % % % ONE-SIDED FOURIER TRANSFORM — 
function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

.

[Adib Muhammad]

Thank you for your answer! It really help

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[Jaimin]

Hi @Adib Muhammad

Begin by creating a complete time series with your data, ensuring any missing hours are filled in. You can use the “interp1” function to fill in the missing values.

You can utilise “fft” function to convert data to frequency domain.

To utilize “plomb” in MATLAB, you might consider a custom implementation, as MATLAB doesn't have a built-in function specifically for “plomb”. However, you can use functions like “pburg” or similar methods to estimate the power spectrum of unevenly sampled data.

Kindly refer following Code Snippet to understand.

% Perform FFT 
N = length(full_data_vector); 
Y = fft(full_data_vector); 

% Using pburg for estimating the power spectrum 
p = 10; % Order of the autoregressive model 
[pxx, f] = pburg(full_data_vector, p, length(full_data_vector), 1); % 1 Hz sampling rate 

Kindly refer following MathWorks Documentations to understand more on

“fft” (https://www.mathworks.com/help/matlab/ref/fft.html),

“pburg” (https://www.mathworks.com/help/signal/ref/pburg.html),

“interp1” (https://www.mathworks.com/help/signal/ref/pburg.html),

“posixtime” (https://www.mathworks.com/help/matlab/ref/datetime.posixtime.html) ,

“datatime” (https://www.mathworks.com/help/matlab/ref/datetime.html),

“readtable” (https://www.mathworks.com/help/matlab/ref/readtable.html).

I hope this will be helpful.

[Adib Muhammad]

thank you for your answer