Is there any way to do LFIT (general linear least-squares fit ) on the time series data using MATLAB?

How to do the attched linear least-square fitting equation on time series data in matlab?

y = a + b*t + c*t^2 + d*sin(2*pi*t) + e*cos(2*pi*t) + f*sin(2*pi*t*2) + g*cos(2*pi*t*2)

a,b,c ... are parameters. t is time and y is dependent variable.

Answer 1

[Star Strider]

That depends entirely on what you want to do.

One approach:

syms a b c d e f g t 
y = a + b*t + c*t^2 + d*sin(2*pi*t) + e*cos(2*pi*t) + f*sin(2*pi*t*2) + g*cos(2*pi*t*2);
yfcn = matlabFunction(y, 'Vars',{[a b c d e f g],t})

yfcn = function_handle with value:
    @(in1,t)in1(:,1)+in1(:,2).*t+in1(:,5).*cos(t.*pi.*2.0)+in1(:,7).*cos(t.*pi.*4.0)+in1(:,4).*sin(t.*pi.*2.0)+in1(:,6).*sin(t.*pi.*4.0)+in1(:,3).*t.^2

This will work in the Optimization Toolbox and Statictics and Machine Learning Toolbox nonlinear parameter estimation funcitons. Specify the initial parameter estimates as a row vector.

tv = linspace(0, 10, 50);                                       % Create Data
yv = exp(-(tv-5).^2) + randn(size(tv))/10;                      % Create Data
B0 = rand(1,7);                                                 % Choose Appropriate Initial Parameter Estimates

B = nlinfit(tv(:), yv(:), yfcn, B0)

B = 1×7

-0.1691 0.2401 -0.0259 -0.0227 -0.0012 0.0266 -0.0076

The Curve Fitting Toolbox does it differently:

yfit = fittype('a + b*t + c*t^2 + d*sin(2*pi*t) + e*cos(2*pi*t) + f*sin(2*pi*t*2) + g*cos(2*pi*t*2)', 'independent',{'t'});
yfit = fit(tv(:) ,yv(:), yfit,'start',B0)

yfit = 
     General model:
     yfit(t) = a + b*t + c*t^2 + d*sin(2*pi*t) + e*cos(2*pi*t) + f*sin(2*pi*t*2)
                     + g*cos(2*pi*t*2)
     Coefficients (with 95% confidence bounds):
       a =     -0.1691  (-0.4125, 0.07436)
       b =      0.2401  (0.1275, 0.3526)
       c =    -0.02595  (-0.03684, -0.01506)
       d =    -0.02275  (-0.1431, 0.09761)
       e =   -0.001153  (-0.1193, 0.117)
       f =     0.02665  (-0.09346, 0.1468)
       g =   -0.007593  (-0.1257, 0.1105)

Experiment to get the result you want.

.

[DIPENDRA007]

@Star Strider Thanks for your help. Is there any optimization technique available to find out the best suited coeffecients of the fitted model ?. The equation doesn't fits with my actual data and giving erroneous values. I have attached the data for your reference.

Thanks in advance.

[Star Strider]

One option is to create as an actual linear regression design matrix —

syms a b c d e f g t y 
y = a + b*t + c*t^2 + d*sin(2*pi*t) + e*cos(2*pi*t) + f*sin(2*pi*t*2) + g*cos(2*pi*t*2)

y = 

[DMs,trms] = coeffs(y,[a b c d e f g])

DMs = 

trms = 

DMf = matlabFunction(DMs)                                       % Design Matrix

DMf = function_handle with value:
    @(t)[1.0,t,t.^2,sin(t.*pi.*2.0),cos(t.*pi.*2.0),sin(t.*pi.*4.0),cos(t.*pi.*4.0)]

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/696424/dataa.txt')

T1 = 4383×2 table
    Var1     Var2 
    ____    ______

      1        NaN
      2        NaN
      3        NaN
      4     375.12
      5     374.13
      6     375.96
      7     377.11
      8     376.48
      9     376.04
     10     375.64
     11     375.23
     12        NaN
     13     375.41
     14     374.14
     15     374.59
     16     374.99

Lv = ~isnan(T1{:,2});
x = T1{Lv,1};
y = T1{Lv,2};
for k = 1:numel(x)
    DM(k,:) = DMf(x(k));                            % Create Design Matrix From Function
end
Q1 = DM

Q1 = 3952×7

1.0000 4.0000 16.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 5.0000 25.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 6.0000 36.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 7.0000 49.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 8.0000 64.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 9.0000 81.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 10.0000 100.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 11.0000 121.0000 -0.0000 1.0000 -0.0000 1.0000 1.0000 13.0000 169.0000 0.0000 1.0000 0.0000 1.0000 1.0000 14.0000 196.0000 -0.0000 1.0000 -0.0000 1.0000

Bv = DM \ y

Warning: Rank deficient, rank = 3, tol = 4.769522e-04.

Bv = 7×1

0 0.0057 -0.0000 0 0 0 370.5335

figure
plot(x, y, '.')
hold on
plot(x, DM*Bv, '-r')
hold off
grid
xlabel('x')
ylabel('y')
title('Linear Approximation')

Fs = 1;                                                     % Sampling Frequency
[yr,xr] = resample(y,x,Fs);                                 % Resample Data
Fn = Fs/2;                                                  % Nyquist Frequency
yrm = y - mean(y);                                          % Remove Mean To Show Peaks
L = numel(xr);                                              % Original Matrix Length
N = 2^nextpow2(L);                                          % FFT Length
FTyr = fft(yrm,N)/L;                                        % FFT
Fv = linspace(0, 1, N/2)*Fn;                                % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector

[pks,locs] = findpeaks(abs(FTyr(Iv))*2, 'MinPeakProminence',5);

figure
plot(Fv, abs(FTyr(Iv))*2)
hold on
plot(Fv(locs), pks, '^r')
hold off
grid
xlim([0 0.01])
text(Fv(locs),pks,compose('$\\leftarrow Frequency\\ %.5f$ \n $\\ \\ \\ \\ Amplitude\\ %.1f$',[Fv(locs).' pks]), 'Horiz','left', 'Vert','top', 'Interpreter','latex')

It will be necessary to add terms to vary the frequencies in the trigonometric terms in the original equation, as well as the existing terms in order to accurately estimate the coefficients. The Fourier transform here will give initial parameter estimates for the frequencies and amplitudes. Then use one of the nonllinear aproaches I outlined in my original Answer, along with appropriate other terms from the linear approximation (‘Bv’) in order to estimate all of them.

I leave that to you.

I will help if you have problems with it.

.

[DIPENDRA007]

@Star Strider thanks I have calculated the coeffecients. One more thing, Can we filter the residuals of the original data about the function?

Thnks in advance.

[Star Strider]

My pleasure.

I am not certain what you intend by ‘filter the residuals’. I would just leave them as they are. They are part of the data.

What coefficients did you estimate, and what objective function did you use in the fit?

.

[DIPENDRA007]

@Star Strider residuals (t)=original data(t)- fitted equation data(t). The residuals are filtered to define short term variations and interannual variations that are not determined by the polynomial terms of the equations as well as for smoothing. Can we filter the residual two times using the attached low pass filter equation?

Where fc is chosen such that h(f) = 0.5 at f=fc.

thanks in advance.

[Star Strider]

I still do not understand what you want to do.

If you want to filter the data, use the lowpass (or bandpass, to eliminate the baseline offset and drift) functions. I have no idea what you want to do, however

[yr_filt,lpdf] = lowpass(yr, 4E-3, Fs, 'ImpulseResponse','iir');

will work for the lowpass filter, and

[yr_filt,bpdf] = bandpass(yr,[1E-5, 4E-3], Fs, 'ImpulseResponse','iir');

for the bandpass design (although it might be necessary to experiment with the lower passband edge).

.