MATLAB Answers

Zero values of a fitted curve

47 views (last 30 days)
ludovico
ludovico on 25 Sep 2015
Commented: ludovico on 5 Oct 2015
Hi to everybody, i have a time history dataset that i have fitted with a fourier series to find its peaks. Now i am in need of the zero crossing values of the curve, and as output the corresponding values of time of these zero crossing. Is there a function that does it automatically? I have tried with fnzeros but seems not to work for fitted curves but only for functions, am i wrong? Maybe i'm missing some statement in the code. Waiting for someone's reply i would thank in advance those willing to help me :D
  5 Comments
Star Strider
Star Strider on 25 Sep 2015
Without actually seeing your data and the analysis you did, I cannot provide a definitive response.

Sign in to comment.

Accepted Answer

arich82
arich82 on 25 Sep 2015
There's a paper by Prof. Boyd from the University of Michigan that appears to address this:
"Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding" http://link.springer.com/article/10.1007%2Fs10665-006-9087-5#page-1
I think this does what you want (not that the roots of the fourier series are given in the variable t):
% fourier coefficients
% (dummy data)
a = [1, 1, 1, 1];
b = [0, 0, 1, 0];
% Note:
% a(1) corresponds to a_0 in the reference
% a(2:end) corresponds to a(1:N) in the reference;
N = numel(a) - 1;
h = [a(end:-1:2) + 1i*b(end:-1:2), 2*a(1), (a(2:end) - 1i*b(2:end))];
B = diag(ones(2*N - 1, 1), 1);
B(end, :) = -h(1:end - 1)/(a(end) - 1i*b(end));
z = eig(B);
t = angle(z);
% note: duplicate t's are imgaginary roots
% eliminate duplicates within tolerance
tol = 10*eps;
t = sort(t);
ind = find(diff(t)./abs(t(2:end)) < tol);
ind = [ind; ind + 1];
t(ind(:)) = [];
% plot results
% f = ugly inline fourier series anonymous function
f = @(x, a, b) sum(bsxfun(@times, a(:).', cos(x(:)*(0:numel(a) - 1))) + bsxfun(@times, b(:).', sin(x(:)*(0:numel(b)-1))), 2);
x = linspace(-pi, pi, 128 + 1);
hf = figure;
ha = axes;
axis([-pi, pi, -sum(a+b)-0.5, sum(a+b)+0.5]);
hold(ha, 'on');
plot(x, f(x, a, b));
plot(t, f(t, a, b), 'o');
plot([-pi, pi], [0, 0], '--c');
I probably won't be around much in the next two days, but feel free to leave a comment, and I'll try to get back to you.

More Answers (2)

Matt J
Matt J on 25 Sep 2015
Edited: Matt J on 25 Sep 2015
fnzeros only appears to work with spline fits, which makes sense. That's the only way the zeros will have an analytical closed form solution. However, one option would be to sample your Fourier series fit at a finely spaced selection of points. Then, you can fit a spline to those points and use fnzeros on the spline fit to approximate the zeros of the Fourier fit.
If needed, you can refine further by using fzero() on the Fourier fit to do a numerical search. The result of fnzeros from above will give you a list of the initial seed points or search intervals that you would loop over, applying fzero() to each one.

ludovico
ludovico on 2 Oct 2015
Guys i must confess i'm completely lost with this issue, i will post here my code in hope that somebody would help me understand what to do. As i already tried to explain, i have a time series of data, i have fitted it with fourier and now, on this fitted curve that is sinusoidal and symmetric over the zero value, i need to find the zero crossing and the exact location of the time of these crossing. Hope to have been clear in the explanation this time! Thanks to everybody by the way!
clear all
close all
clc
% INPUT per fitting Lift%
[TH] = load('Cl_OR2_f0a_4_N.txt');
% assign first column to time step and second to Cl%
Time_step = TH(:,1);
Cl = TH(:,2);
%calc pos e neg peaks%
[pks_1,locs_1] = findpeaks(Cl,'MINPEAKDISTANCE',20);
[pks_2,locs_2] = findpeaks(-Cl,'MINPEAKDISTANCE',20);
%plot function Cl with blu%
f1 = figure;
plot(Time_step,Cl,'b');
%set dim grid assis x to 1sec%
set(gca,'xtick',[20:1:40]);
set(gca,'Ytick',[-2.5:0.1:2.5]);
hold on;
grid on;
%plot function peaks in red%
plot(Time_step(locs_1),pks_1,'ro');
plot(Time_step(locs_2),-pks_2,'ro');
grid on;
%calc mean pos peaks%
M1 = mean(pks_1)
%calc val(RMS)%
Cl_RMS=0.707*M1
%fitting with Fourier series to derive parameters a0 - a1- b1 %
[fitobject]=fit(Time_step,Cl,'fourier1')
a0=fitobject.a0;
a1=fitobject.a1;
b1=fitobject.b1;
w=fitobject.w;
%calc phase signal%
PhaseCl = atan(a1/b1)
%open curve fitting tool for Fourier coeff calculation with GUI - deactivated%
%cftool(Time_step,Cl)%
Cl_v1 = Cl_RMS*sin(PhaseCl)
Cl_a1 = Cl_RMS*(-cos(PhaseCl))
%declaration functions SPOST - VEL - ACC for the set freq%
A= 0.015;
f=0.53;
t=[20:0.1:40];
SPOST = (A/2)*sin(2*pi*f*Time_step);
VEL = (A/2)*(2*pi*f)*cos(2*pi*f*Time_step);
ACC = (-A/2)*((2*pi*f)^2)*sin(2*pi*f*Time_step);
plot(Time_step,SPOST,'b');
hold on
plot(Time_step,VEL,'g');
hold on
plot(Time_step,ACC,'color',[0 0.5 0]);
legend('Displacement (m)','Velocity (m/s)','Acceleration (m/sec^2)');
grid on
[fitobject2]=fit(Time_step,SPOST,'fourier1')
a0=fitobject2.a0;
a1=fitobject2.a1;
b1=fitobject2.b1;
w=fitobject2.w;
PhaseSPOST = atan(a1/b1)
Signal_phase = PhaseCl - PhaseSPOST
  3 Comments
ludovico
ludovico on 5 Oct 2015
Solved the issue mate, many many thanks really!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by