convert time domain data to frequency domain
Afficher commentaires plus anciens
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!
1 commentaire
Adib Muhammad
le 29 Oct 2024
Réponse acceptée
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'}
% 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)
dt = diff(hour(T1{:,1}));
nonseq = nnz(dt>1) % Number Of Non-Sequantial Times
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))
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))
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 commentaires
Adib Muhammad
le 30 Oct 2024
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'}
% 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)
dt = diff(hour(T1{:,1}));
nonseq = nnz(dt>1) % Number Of Non-Sequantial Times
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))
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))
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
le 1 Nov 2024
Star Strider
le 1 Nov 2024
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
Plus de réponses (1)
Jaimin
le 29 Oct 2024
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
I hope this will be helpful.
1 commentaire
Adib Muhammad
le 30 Oct 2024
Catégories
En savoir plus sur Filter Analysis dans Centre d'aide et File Exchange
le 29 Oct 2024
le 1 Nov 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
