How to estimate annual and semi-annual signals?

8 vues (au cours des 30 derniers jours)
hanif hamden
hanif hamden le 25 Nov 2020
Commenté : Mathieu NOE le 26 Nov 2020
Hi everyone, I would like to estimate annual ans semi-annual signal from my data. Attached is my data as well as my coding, I managed to estimate bias and linear trend. However, I was stuck in estimating annual and semi-annual signal. I hope anyone can help me on this matter. Thank you.
A = load('FPT1.txt');
t = A(:,1);
slv = A(:,5);
dt1 = num2str(t,'%.2f');
dtm1 = datetime(dt1,'InputFormat','yyyyMMddHHmmss.SS');
dtn1 = datenum(dtm1);
dtmFill1 = fillmissing(dtm1,'linear');
%%PARAMETER FIT
% hobs = h0 + h1T + h2*sin(2*pi*T) + h3*cos(2*pi*T) + h4*sin(4*pi*T) + h5*cos(4*pi*T)
% Estimate bias/offset
h0 = mean(slv);
% Estimate linear trend
b = polyfit(dtn1,slv,1);
h1 = polyval(b, dtn1);
%% Stuck here
% Estimate Annual Signals [h2*sin(2*pi*T) + h3*cos(2*pi*T)]
% However, the time for this data is every 9.9156 days
% f1=1/(365*24); %frequency of 1 yr oscillations
h2 = sin(2*pi*T)/slv;
h3 = cos(2*pi*T)/slv;
% Estimate Semi-Annual Signals [h4*sin(4*pi*T) + h5*cos(4*pi*T)]
% f2=1/((365*24))/2; %frequency of 1/2 year oscillations
h4 = slv*sin(4*pi*T);
h5 = slv*cos(4*pi*T);

Réponse acceptée

Mathieu NOE
Mathieu NOE le 25 Nov 2020
hello
you're not stuck anymore
below my code (the h parameters computation is the same as DFT coefficients)
if you're stisfied, pls accept my answer
all the best
A = load('FPT1.txt');
t = A(:,1);
slv = A(:,5);
dt1 = num2str(t,'%.2f');
dtm1 = datetime(dt1,'InputFormat','yyyyMMddHHmmss.SS');
dtn1 = datenum(dtm1);
% resample with constant time sampling
% remember time is in days
dt = 7; % 7 days interval
Fs = 1/dt; % NB not in Hz = 1/s but in 1/day
new_t = min(dtn1):dt:max(dtn1);
new_data = interp1(dtn1,slv,new_t);
% figure(1),plot(dtn1,slv,'b+-',new_t,new_data,'r*-');
%%PARAMETER FIT
% hobs = h0 + h1T + h2*sin(2*pi*T) + h3*cos(2*pi*T) + h4*sin(4*pi*T) + h5*cos(4*pi*T)
% Estimate bias/offset
% h0 = mean(new_data);
% Estimate linear trend
b = polyfit(new_t,new_data,1);
% h1 = polyval(b, new_t);
%% Not anymore Stuck here
% Estimate Annual Signals [h2*sin(2*pi*T) + h3*cos(2*pi*T)]
% However, the time for this data is every 9.9156 days
f1=1/(365); %frequency of 1 yr oscillations time base in days)
h2 = 2/length(new_t)*sum(sin(2*pi*f1*new_t).*new_data);
h3 = 2/length(new_t)*sum(cos(2*pi*f1*new_t).*new_data);
% Estimate Semi-Annual Signals [h4*sin(4*pi*T) + h5*cos(4*pi*T)]
f2=2*f1; %frequency of 1/2 year oscillations
h4 = 2/length(new_t)*sum(sin(2*pi*f2*new_t).*new_data);
h5 = 2/length(new_t)*sum(cos(2*pi*f2*new_t).*new_data);
% reconstructed signal
hobs = h2*sin(2*pi*f1*new_t) + h3*cos(2*pi*f1*new_t) + h4*sin(2*pi*f2*new_t) + h5*cos(2*pi*f2*new_t);
% use the delta = new_data-hobs, for better linear trend estimation
delta = new_data-hobs;
% Estimate linear trend
b = polyfit(new_t,delta,1);
h0 = b(2);
h1 = b(1);
hobs = h0 + h1*new_t + h2*sin(2*pi*f1*new_t) + h3*cos(2*pi*f1*new_t) + h4*sin(2*pi*f2*new_t) + h5*cos(2*pi*f2*new_t);
figure(2),plot(dtn1,slv,'b+-',new_t,hobs,'r-*');
  2 commentaires
hanif hamden
hanif hamden le 26 Nov 2020
Thank you sir. the code works well. Just one more question. I wanted to subtract the hobs model with the original data, since the matrix is not the same. Is it correct if I interpolate the hobs based on time?
% -----Additional---
new_t2 = new_t';
hobs2 = hobs';
dtm2 = datetime(new_t2,'ConvertFrom','datenum');
Ts_hobs = [new_t2 hobs2];
new_Ts_hobs = interp1(Ts_hobs(:,1),Ts_hobs(:,2),dtn1,'makima');
Nslv = slv-new_Ts_hobs;
figure(3),plot(dtn1,slv,'b+-',dtn1,Nslv,'r-*');
Mathieu NOE
Mathieu NOE le 26 Nov 2020
yes it's ok
just simplified a bit your code;
this is doing the same with less :
% -----Additional---
dtm2 = datetime(new_t,'ConvertFrom','datenum');
new_Ts_hobs = interp1(new_t,hobs,dtn1,'makima');
Nslv = slv-new_Ts_hobs;
figure(3),plot(dtn1,slv,'b+-',dtn1,Nslv,'r-*');

Connectez-vous pour commenter.

Plus de réponses (0)

Produits


Version

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by