Difference between revisions of "Fourier fit"
From Rizzo_Lab
| (One intermediate revision by one other user not shown) | |||
| Line 1: | Line 1: | ||
function [Yp, Parm] = fourier_fit(Ye,X,N) | function [Yp, Parm] = fourier_fit(Ye,X,N) | ||
% written by Trent E. Balius | % written by Trent E. Balius | ||
| + | % modified by Lingling Jiang | ||
% Robert C. Rizzo Research Group | % Robert C. Rizzo Research Group | ||
% Stony Brook University | % Stony Brook University | ||
| Line 8: | Line 9: | ||
% Parm = {V,A1,B1,A2,B2} | % Parm = {V,A1,B1,A2,B2} | ||
% f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)} | % f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)} | ||
| + | % f(x) = P(1) + sum(i=1:N){P(2*i)*cos(i*x)+P(2*i+1)*sin(i*x)} | ||
| + | |||
| + | % f(x) = P(1) + sum(j=2:2:2*N){P(j)*cos((j/2)*x)} | ||
| + | % +sum(j=3:2:2*N+1){P(j)*sin(((j-1)/2)*x)} | ||
| + | |||
% Ye is the Y value from the data | % Ye is the Y value from the data | ||
% X is the X value from the data | % X is the X value from the data | ||
| Line 30: | Line 36: | ||
%whos | %whos | ||
% f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)} | % f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)} | ||
| + | % f(x) = P(1) + sum(i=1:N){P(2*i)*cos(i*x)+P(2*i+1)*sin(i*x)} | ||
| + | |||
| + | % f(x) = P(1) + sum(j=2:2:2*N){P(j)*cos((j/2)*x)} | ||
| + | % +sum(j=3:2:2*N+1){P(j)*sin(((j-1)/2)*x)} | ||
Yp = Parm(1)*ones(size(X)); | Yp = Parm(1)*ones(size(X)); | ||
| Line 35: | Line 45: | ||
for i = 2:length(Parm) | for i = 2:length(Parm) | ||
if (mod(i,2) == 0) % i is even | if (mod(i,2) == 0) % i is even | ||
| − | Yp = Yp + Parm(i)*cos(i*X); | + | Yp = Yp + Parm(i)*cos((i/2)*X); |
else % i is odd | else % i is odd | ||
| − | Yp = Yp + Parm(i)*sin(i*X); | + | Yp = Yp + Parm(i)*sin((i-1)/2*X); |
end | end | ||
end | end | ||
return | return | ||
Latest revision as of 07:23, 6 May 2011
function [Yp, Parm] = fourier_fit(Ye,X,N)
% written by Trent E. Balius
% modified by Lingling Jiang
% Robert C. Rizzo Research Group
% Stony Brook University
% syntax: [Yp,Parm] = fourier_fit(Ye,X,N)
% Yp is the predicted Y values
% Parm is the parameter set = 2*N + 1
% Parm = {V,A1,B1,A2,B2}
% f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)}
% f(x) = P(1) + sum(i=1:N){P(2*i)*cos(i*x)+P(2*i+1)*sin(i*x)}
% f(x) = P(1) + sum(j=2:2:2*N){P(j)*cos((j/2)*x)}
% +sum(j=3:2:2*N+1){P(j)*sin(((j-1)/2)*x)}
% Ye is the Y value from the data
% X is the X value from the data
% N is the number over which sum is calculated
% uses fmincon
input = rand( (2*N + 1),1); %intial guess
lb = -1000*ones(2*N + 1,1); % lower bound
ub = 1000*ones(2*N + 1,1); % upper bound
[Parm,f] = fmincon(@objective_function,input,[],[],[],[],lb,ub,[],[],Ye,X);
Yp = fourier_series(Parm,X);
function [val] = objective_function(Parm,Y,X)
Yp = fourier_series(Parm,X);
%[Y,Yp]
val = sum((Y-Yp).^2);
return
function [Yp] = fourier_series(Parm,X)
%whos
% f(x) = V(1) + sum(i=1:N){A1(i)*cos(i*x)+B1(i)*sin(i*x)}
% f(x) = P(1) + sum(i=1:N){P(2*i)*cos(i*x)+P(2*i+1)*sin(i*x)}
% f(x) = P(1) + sum(j=2:2:2*N){P(j)*cos((j/2)*x)}
% +sum(j=3:2:2*N+1){P(j)*sin(((j-1)/2)*x)}
Yp = Parm(1)*ones(size(X));
for i = 2:length(Parm)
if (mod(i,2) == 0) % i is even
Yp = Yp + Parm(i)*cos((i/2)*X);
else % i is odd
Yp = Yp + Parm(i)*sin((i-1)/2*X);
end
end
return
