フィルターのクリア

How to estimate annual and semi-annual signals?

4 ビュー (過去 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 ExchangeSmoothing and Denoising についてさらに検索

製品


リリース

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by