1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Mar 5, 2013 #1
    1. The problem statement, all variables and given/known data

    screen-capture-3-37_zps35a0465b.png

    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 (?).

    fig1_zps4d964c93.jpg

    fig2_zps8eb40594.jpg


    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))

    2. Relevant 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

    3. 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.
     
  2. jcsd
  3. Mar 5, 2013 #2

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    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?
     
  4. Mar 5, 2013 #3
    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={c0, c1, ..., c30}, and I can evaluate it at the larger set of data points using polyval like I did.


    % 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;
     
  5. Mar 6, 2013 #4

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    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.

    fsum is a function? Oh, okay.
    I don't understand why you need two different cases here.
     
  6. Mar 6, 2013 #5
    Could you explain the second bullet point to me?
     
  7. Mar 7, 2013 #6

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    Why do you use different formulas for k<=15 and k>15?
     
  8. Mar 7, 2013 #7
    Isn't that what it shows on the assignment?
     
  9. Mar 8, 2013 #8

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    Right, but I don't understand the reason there either. It might be related to the specific way the fourier transform is calculated.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Fitting a curve using a spline, Fourier transform, etc.
Loading...