Fitting a curve using a spline, Fourier transform, etc.

[Jamin2112]

Homework Statement

Just wondering if my output seems wrong. The interpolating polynomial looks like it's way off, though I've looked over my code many times and it seems right (?).

clc
clear all
format long

x1=[1:1/10:4];
y1=zeros(1,length(x1));
y1(find(x1<=3))=floor(x1((find(x1<=3))).^2);
y1(find(x1>3))=9*cos(3*x1(find(x1>3))-9).^3;
xx=[1:1/1000:4];
yy=zeros(1,length(xx));
yy(find(xx<=3))=floor((xx(find(xx<=3))).^2);
yy(find(xx>3))=9*cos(3*xx(find(xx>3))-9).^3;
p=polyfit(x1,y1,length(x1)+1);
y30=polyval(p,xx);
ffty1=fft(y1);
fsum=0;
for k=0:15
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*k*(xx-1)/3);
end
for k=16:30
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*(k-31)*(xx-1)/3);
end
f31=fsum/31;
ynakstruct=spline(xx,yy);
ynakcoeffs=ynakstruct.coefs;
ynak=zeros(1,length(xx));
ynak(1)=yy(1);
for k=2:length(ynak)
ynak(k)=ynakcoeffs(k-1,4)+ynakcoeffs(k-1,3)*(xx(k)-xx(k-1))+ynakcoeffs(k-1,2)*(xx(k)-xx(k-1))^2+ynakcoeffs(k-1,1)*(xx(k)-xx(k-1))^3;
end
yclstruct=spline(xx,[0 yy 0]);
yclcoeffs=yclstruct.coefs;
ycl=zeros(1,length(xx));
ycl(1)=yy(1);
for i=2:length(ycl)
ycl(k)=yclcoeffs(k-1,4)+yclcoeffs(k-1,3)*(xx(k)-xx(k-1))+yclcoeffs(k-1,2)*(xx(k)-xx(k-1))^2+yclcoeffs(k-1,1)*(xx(k)-xx(k-1))^3;
end
figure;
subplot(2,2,1)
hold on
plot(xx,yy,xx,y30,'-')
legend('yy', 'y30')
e1=(1/3001)*norm(y30-yy,1);
title(sprintf('Error: %d', e1))
subplot(2,2,2)
hold on
plot(xx,yy,xx,f31,'-')
legend('yy', 'f31')
e2=(1/3001)*norm(f31-yy,1);
title(sprintf('Error: %d', e2))
subplot(2,2,3)
hold on
plot(xx,yy,xx,ynak,'-')
legend('yy', 'ynak')
e3=(1/3001)*norm(ynak-yy,1);
title(sprintf('Error: %d', e3))
subplot(2,2,4)
hold on
plot(xx,yy,xx,ycl,'-')
legend('yy', 'ycl')
e4=(1/3001)*norm(ycl-yy,1);
title(sprintf('Error: %d', e4))

x2=[1:1/100:4];
y2=zeros(1,length(x2));
y2(find(x2<=3))=floor(x2((find(x2<=3))).^2);
y2(find(x2>3))=9*cos(3*x2(find(x2>3))-9).^3-9;
yy=zeros(1,length(xx));
yy(find(xx<=3))=floor((xx(find(xx<=3))).^2);
yy(find(xx>3))=9*cos(3*xx(find(xx>3))-9).^3;
p=polyfit(x2,y2,length(x2)+1);
y301=polyval(p,xx);
ffty2=fft(y2);
fsum=0;
for k=0:150
fsum=fsum+(ffty2(k+1))*(k+1)*exp(2*pi*i*k*(xx-1)/3);
end
for k=151:300
fsum=fsum+(ffty2(k+1))*(k+1)*exp(2*pi*i*(k-301)*(xx-1)/3);
end
f301=fsum/301;
ynakstruct=spline(xx,yy);
ynakcoeffs=ynakstruct.coefs;
ynak=zeros(1,length(xx));
ynak(1)=yy(1);
for k=2:length(ynak)
ynak(k)=ynakcoeffs(k-1,4)+ynakcoeffs(k-1,3)*(xx(k)-xx(k-1))+ynakcoeffs(k-1,2)*(xx(k)-xx(k-1))^2+ynakcoeffs(k-1,1)*(xx(k)-xx(k-1))^3;
end
yclstruct=spline(xx,[0 yy 0]);
yclcoeffs=yclstruct.coefs;
ycl=zeros(1,length(xx));
ycl(1)=yy(1);
for k=2:length(ycl)
ycl(k)=yclcoeffs(k-1,4)+yclcoeffs(k-1,3)*(xx(k)-xx(k-1))+yclcoeffs(k-1,2)*(xx(k)-xx(k-1))^2+yclcoeffs(k-1,1)*(xx(k)-xx(k-1))^3;
end
figure(2)
subplot(2,2,1)
hold on
plot(xx,yy,xx,y301,'-')
legend('yy', 'y301')
e1=(1/3001)*norm(y301-yy,1);
title(sprintf('Error: %d', e1))
subplot(2,2,2)
hold on
plot(xx,yy,xx,f301,'-')
legend('yy', 'f301')
e2=(1/3001)*norm(f301-yy,1);
title(sprintf('Error: %d', e2))
subplot(2,2,3)
hold on
plot(xx,yy,xx,ynak,'-')
legend('yy', 'ynak')
e3=(1/3001)*norm(ynak-yy,1);
title(sprintf('Error: %d', e3))
subplot(2,2,4)
hold on
plot(xx,yy,xx,ycl,'-')
legend('yy', 'ycl')
e4=(1/3001)*norm(ycl-yy,1);
title(sprintf('Error: %d', e4))

Homework Equations

fft(x) returns the Fourier coefficients for the function x given by sample points
polyfit(x,y,n) returns n polynomial coefficients for a polynomial fitted to the points given by the vectors x and y
polyval(p,x) returns the values for a polynomial with coefficients given in the vector p at the points in the vector x

The Attempt at a Solution

So my spline things seem to be working out. It's just the other two that are looking way off, even though I've checked my code a bunch. Just wondering if maybe you can think of some numerical reason why this might be the case.

 

[mfb]

Can you get some debug data from the poly fit?
Did you test it with an easier function? Or with more sampling points? The "fit" has as many equations as the polynomial has degrees of freedom, so it is not a real "fit" here.

What happens to ffty1, and where do you plot it?

 

[Jamin2112]

Can you get some debug data from the poly fit?
Did you test it with an easier function? Or with more sampling points? The "fit" has as many equations as the polynomial has degrees of freedom, so it is not a real "fit" here.

Isn't it asking for a polynomial of order 30? So I should be getting those coefficients from the smaller set of 31 data points y1, for a polynomial of order 30 with coefficients p={c₀, c₁, ..., c₃₀}, and I can evaluate it at the larger set of data points using polyval like I did.

What happens to ffty1, and where do you plot it?

% Get the Fourier coefficients (complex #'s) for the 31 pts of data
ffty1=fft(y1);
% Use those to construct the function f31 at the pts xx={1, 1.01, 1.02, ... 4.0}
fsum=0;
for k=0:15
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*k*(xx-1)/3);
end
for k=16:30
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*(k-31)*(xx-1)/3);
end
f31=fsum/31;

 

[mfb]

Isn't it asking for a polynomial of order 30? So I should be getting those coefficients from the smaller set of 31 data points y1, for a polynomial of order 30 with coefficients p={c₀, c₁, ..., c₃₀}, and I can evaluate it at the larger set of data points using polyval like I did.

Sure, but you use 31 data points to determine 31 coefficients. That is just a linear equation system, and it should have a unique solution, where the polynomial hits all 31 points.

% Get the Fourier coefficients (complex #'s) for the 31 pts of data
[...]
f31=fsum/31;

fsum is a function? Oh, okay.
I don't understand why you need two different cases here.

 

[Jamin2112]

fsum is a function? Oh, okay.
I don't understand why you need two different cases here.

Could you explain the second bullet point to me?

 

[mfb]

for k=0:15
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*k*(xx-1)/3);
end
for k=16:30
fsum=fsum+(ffty1(k+1))*(k+1)*exp(2*pi*i*(k-31)*(xx-1)/3);
end

Why do you use different formulas for k<=15 and k>15?

 

[Jamin2112]

Why do you use different formulas for k<=15 and k>15?

Isn't that what it shows on the assignment?

 

[mfb]

Right, but I don't understand the reason there either. It might be related to the specific way the Fourier transform is calculated.