0001 function C = fourier_fit(points,N,start);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 if ischar(points) && strcmp(points,'UNIT_TEST'); do_unit_test; return ; end
0020
0021 if size(points,2)==2
0022 if nargin<2; N= size(points,1); end
0023 C = fit_to_fourier(points,N);
0024 elseif size(points,2)==1
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
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
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));
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;
0044
0045 function xy = fit_from_fourier(C,linear_frac,start);
0046
0047 N = length(C);
0048 n2 = ceil(N/2);
0049
0050 pad = zeros(10000,1);
0051 if rem(N,2)==0
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
0059
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
0069 zos(end+1) = zos(1);
0070 dpath= abs(diff(zos));
0071 pathlen = [0;cumsum(dpath)];
0072
0073
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);