fourier_fit

PURPOSE ^

FOURIER_FIT: use fourier series to interpolate onto a boundary

SYNOPSIS ^

function C = fourier_fit(points,N,start);

DESCRIPTION ^

 FOURIER_FIT: use fourier series to interpolate onto a boundary

 [pp] = fourier_fit(points) fits a Fourier series
    points - [x y] contour to be fitted
 [pp] = fourier_fit(points,N) fits a Fourier series and downsamples
    N is the number of Fourier components to fit to.

 [xy] = fourier_fit(pp,  linear_frac, start) returns points at linear_frac
 distance along the contour
    pp          - piecewise polynomial structure
    linear_frac - vector of length fractions (0 to 1) to calculate points
    start       - (optional) a seed for the first point
    xy          - cartesian coordinates of the points

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function C = fourier_fit(points,N,start);
0002 % FOURIER_FIT: use fourier series to interpolate onto a boundary
0003 %
0004 % [pp] = fourier_fit(points) fits a Fourier series
0005 %    points - [x y] contour to be fitted
0006 % [pp] = fourier_fit(points,N) fits a Fourier series and downsamples
0007 %    N is the number of Fourier components to fit to.
0008 %
0009 % [xy] = fourier_fit(pp,  linear_frac, start) returns points at linear_frac
0010 % distance along the contour
0011 %    pp          - piecewise polynomial structure
0012 %    linear_frac - vector of length fractions (0 to 1) to calculate points
0013 %    start       - (optional) a seed for the first point
0014 %    xy          - cartesian coordinates of the points
0015 
0016 % (C) Andy Adler, 2010. Licenced under GPL v2 or v3
0017 % $Id: fourier_fit.m 4927 2015-05-07 23:10:04Z aadler $
0018 
0019 if isstr(points) && strcmp(points,'UNIT_TEST'); do_unit_test; return ; end
0020 
0021 if size(points,2)==2 % calling analysis function
0022    if nargin<2; N= size(points,1); end
0023    C = fit_to_fourier(points,N);
0024 elseif size(points,2)==1 % calling synthesis function
0025    if nargin<3; start = []; end
0026    C = fit_from_fourier(points,N,start);
0027 else
0028    error('size of first parameter to fourier_fit not undersood');
0029 end
0030 
0031 % this will crash if N>length(points)
0032 function C = fit_to_fourier(points,N);
0033    z = points*[1;1i];
0034    Z = fft(z,max(N,size(points,1)));
0035    if rem(N,2)==0 % Even
0036      N2 = N/2;
0037      Zlp = Z([1,2:N2,1,end-(N2-2:-1:0)]);
0038      Zlp(N2+1) = 0.5*(Z(N2+1) + Z(end-N2+1)); %centre point
0039    else 
0040      N2 = floor(N/2);
0041      Zlp = Z([1,2:N2+1,end-(N2-1:-1:0)]);
0042    end
0043    C = length(Zlp)/length(Z)*Zlp; % Amplitude scaling
0044     
0045 function xy = fit_from_fourier(C,linear_frac,start);
0046    % Step 1: oversample
0047    N = length(C);
0048    n2 = ceil(N/2);
0049 
0050    pad = zeros(10000,1);
0051    if rem(N,2)==0 % even
0052       Zos = [C(1:n2); C(n2+1)/2; pad; C(n2+1)/2; C(n2+2:end)];
0053    else
0054       Zos = [C(1:n2); pad; C(n2+1:end)];
0055    end
0056    Zos = length(Zos)/length(C)*Zos;
0057    zos = ifft(Zos);
0058    % rearrange the points such that they start as close as possible to the
0059    % seed
0060    if ~isempty(start)
0061        if size(start,2) == 1; start = start'; end
0062        start = start*[1; 1i];
0063        dist = abs(zos-start);
0064        [jnk,p] = min(dist);
0065        zos = circshift(zos,-p+1);
0066    end
0067    
0068    % Step 2:
0069    zos(end+1) = zos(1); % make sure the loop is closed
0070    dpath= abs(diff(zos));
0071    pathlen = [0;cumsum(dpath)];
0072 
0073    % interpolate points onto path
0074    npath = pathlen/max(pathlen);
0075    linear_frac= mod(linear_frac,1);
0076    z_int = interp1(npath, zos, linear_frac);
0077 
0078    xy= [real(z_int(:)), imag(z_int(:))];
0079    
0080 function do_unit_test
0081 
0082     a = [
0083    -0.8981   -0.7492   -0.2146    0.3162    0.7935    0.9615    0.6751    0.0565   -0.3635   -0.9745
0084     0.1404    0.5146    0.3504    0.5069    0.2702   -0.2339   -0.8677   -0.6997   -0.8563   -0.4668 ]';
0085 
0086    C= fourier_fit(a,10);
0087    xy = fourier_fit(C, linspace(.05,1.04,100));
0088    xy2= fourier_fit(C, [.3,.4,.5,.6]);
0089    plot(a(:,1),a(:,2),'r',xy(:,1),xy(:,2),'b',xy2(:,1),xy2(:,2),'m*');
0090 
0091    a(5,:)= [];
0092    C= fourier_fit(a);
0093    xy = fourier_fit(C, linspace(.05,1.04,100));
0094    xy2= fourier_fit(C, [.3,.4,.5,.6]);
0095    plot(a(:,1),a(:,2),'r',xy(:,1),xy(:,2),'b',xy2(:,1),xy2(:,2),'m*');
0096    eidors_msg('VIEW GRAPH TO VERIFY',0);

Generated on Tue 12-May-2015 16:00:42 by m2html © 2005