Specifically:
D = load('Louise data_TRT_HELP.mat');
data = D.data_TRT;
t= linspace(0, numel(data), numel(data));
y=data;
Fs=500;
Fn = Fs/2; % Nyquist Frequency
L = size(t,2);
y = detrend(y); % Remove Linear Trend
fty = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
Iv = 1:length(Fv); % Index Vector
figure(1)
plot(Fv, abs(fty(Iv)))
grid
Wp = 30/Fn; % Passband (Normalised)
Ws = 40/Fn; % Stopband (Normalised)
Rp = 1; % Passband Ripple
Rs = 50; % Stopband Ripple
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp); % Transfer Function Coefficients
[sos,g] = zp2sos(z,p,k); % Second-Order-Section For Stability
figure(2)
freqz(sos, 2^14, Fs)
yf = filtfilt(sos,g,y);
yu = max(yf);
yl = min(yf);
yr = (yu-yl); % Range of ‘y’
yz = yf-yu+(yr/2);
zt = t(yz(:) .* circshift(yz(:),[1 0]) <= 0); % Find zero-crossings
per = 2*mean(diff(zt)); % Estimate period
ym = mean(y); % Estimate offset
fit = @(b,x) b(1) .* exp(b(2).*x) .* (sin(2*pi*x./(b(3)+b(4).*x) + 2*pi/b(5))) + b(6); % Objective Function to fit
fcn = @(b) sum((fit(b,t) - yf).^2); % Least-Squares cost function
s = fminsearch(fcn, [yr; -1; per; 0.01; -1; ym]) % Minimise Least-Squares
xp = linspace(min(t),max(t));
figure(3)
plot(t,y,'b')
hold on
plot(t,yf,'g', 'LineWidth',2)
plot(xp,fit(s,xp), 'r', 'LineWidth',1.5)
hold off
grid
title('Sample Data 2')
xlabel('Time')
ylabel('Amplitude')
legend('Original Data', 'Lowpass-Filtered Data', 'Fitted Curve')
producing:
