How to estimate annual and semi-annual signals?

5 次查看（过去 30 天）
hanif hamden
hanif hamden 2020-11-25
评论： Mathieu NOE 2020-11-26
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);

请先登录，再回答此问题。

采纳的回答

Mathieu NOE
Mathieu NOE 2020-11-25
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 个评论
hanif hamden
hanif hamden 2020-11-26
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 2020-11-26
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-*');

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Model Identification 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by