convert time domain data to frequency domain

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

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2024-10-30，9:17

在 MATLAB Online 中打开

Time_domain_data.xlsx

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 2024-11-1，2:17

Thank you for your answer! It really help

Star Strider 2024-11-1，2:28

My pleasure!

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

.

请先登录，再进行评论。

Answer 2

Jaimin 2024-10-29，6:32

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2162650-convert-time-domain-data-to-frequency-domain#answer_1538610

在 MATLAB Online 中打开

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 