Zero values of a fitted curve

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

[Star Strider]

‘...the zero crossing values of the curve...’

What curve? The original time series, the Fourier transform, or something else? How did you express the Fourier transform (for instance, amplitude, or separate real and imaginary components)?

[ludovico]

Firstly thanks for the reply, i need to find the zeros of the fitted curve, represented in harmonic components. Shall i create a function from the fitted curve and than use the fnzeros??

[Star Strider]

I still don’t understand what the ‘fitted curve’ is. If you have an objective function and a set of estimated parameters (such as those estimated by nlinfit or lsqcurvefit), you can use those parameters and the objective function to find the zero-crossings.

[ludovico]

you're right i have a plotted curve of the values of the time history, and i made a fitting with Fourier series to find his peaks. can i find the zeros from the plotted graph or shall i make a function out of the input data?

[Star Strider]

Without actually seeing your data and the analysis you did, I cannot provide a definitive response.

Answer 1

[arich82]

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

A summary seems to be given here: http://math.stackexchange.com/questions/370996/roots-of-a-finite-fourier-series

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.

Answer 2

[Matt J]

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.

Answer 3

[ludovico]

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

[arich82]

It's still not entirely clear what you want:

Apparently, you have analytical expressions for SPOST, VEL, and ACC, and you want to find the zero crossing for, what exactly? (You also declared a t=[20:0.1:40], but never used it...)

I see the fitObject for C1, and another fitObject2 for SPOST. You could use the code from my previous post by defining a=[a0, a1] and b=[0, b1], and the result t would have to be extended periodically over n*2*pi to get the answer in the desired interval, and presumably, rescaled using w. (I don't have the appropriate toolbox, so I don't know the details of the fit function; the online documentation shows a p = 2*pi/(max(xdata)-min(xdata)), but not a w.)

However, for such a simple expression (i.e. a single sin & cos term), you can just compute the result analytically: assuming a0==0 (is this what you mean by "symmetric over the zero value"?), you can combine the sin and cos terms into a single sin (as you appear to have done), and solve for the time t when the entire argument becomes any integer multiple of pi.

For example, if your analytical expression is

 y = 0 + a1*cos(w*t) + b1*sin(w*t)

let,

 A = sqrt(a1^2 + b1^2);
 phi = atan(a1/b1);
 y == A*sin(w*t + phi)

then your roots would be

 t = (pi*[-n:1:n] - phi)/w

for all integer n.

If you're still not clear, you'll have to explain exactly what it is you're having difficulty computing. (Also, it's better to post additional info in the comments, rather than posting it as a new answer...)

[ludovico]

Firstly thanks so much for your availability to my issue, i have to match the displacement, velocity and acceleration of a body with a time series loaded in matlab with this statement following:

    % INPUT per fitting Lift%
    [TH] = load('Cl_OR2_f0a_4_N.txt');

to do so i have to find the fitting of the time series, which is a sinusoidal curve symmetrical on the zero value of the Y assis, and find the crossing of THIS fitted curve. i will link a picture here to explain better the situation. The legend is wrong too, i have to fix it but the green curves are the motion ones, whereas the blue one is the fitted curve i have to find peaks of.

thanks to your last suggestion i've understood something more, i am very grateful to you, i will try it and eventually close this discussion. Sincerely

[ludovico]

Solved the issue mate, many many thanks really!