正交多项式程序共享

土木年华

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('输入文件名:','s');
fid=fopen(fni,'r');
mn=fscanf(fid,'%d',1);
df=fscanf(fid,'%f',1);
fno=fscanf(fid,'%s',1);
h=fscanf(fid,'%f',[2,inf]);
status=fclose(fid);
nm=mn*2;
n=length(h(1,:));
f=0:df:(n-1)*df;
w=[-f(1,n:-1:1) f (1:n)]/max(f);
H=[h(1,n:-1:1) -i*h(2,n:-1:1) h(1,1:n)+i*h(2,1:n)];
[k,1]=size(w);
if k<1
    w=w.';
end
[k,1]=size(H);
if k<1
    H=H.';
end
[p,cp]=fun851(H,w,1,nm);
[q,cq]=fun851(H,w,2,nm);
P=p;
for k=1:nm
    Q(:,k)=(H.*q(:,k));
end
T=-real(P'*Q);
g=real(P'*(H.*q(:,nm+1)));
D=-inv(eye(nm)-T'*T)*T'*g;
C=g-T*D;
D(nm+1)=1;
A=cp*C;
A=A(nm+1:-1:1).';
B=cq*D;
B=B(nm+1:-1:1).';
[R,U,K]=residue(A,B);
F1=abs(U)*max(f);
D1=-real(U)./(abs(U));
gor k=1:nm
if k==1
    u(1:nm-1)=U(2:nm);
else
    u(1:k-1)=U(1:k-1);
    u(k:nm-1)=U(k+1:nm);
end
y=poly(u);
S1(k)=polyval(A,U(k))/polyval(y,U(k));
end
H=h(1,1:n)+i*h(2,1:n);
for l=1:n
    H1(l)=sum(C'.*p(n+1,:))/sum(D'.*q(n+1,:));
end
subplot(2,1,1);
plot(f,real(H),':',f,real(H1),'r');
xlabel('频率(HZ)');
ylabel('实部');
legend('实测''拟合');
grid on;
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1),'r');
xlabel('频率(HZ)');
ylabel('虚部');
legend('实测''拟合');
grid on;
[F2,I]=sort(F1);
m=0;
for k=1:n-1
    if F2(k)~=F2(k+1)
        continue;
    end
    m=m+1;
    l=I(k);
    F(m)=F1(l);
    D(m)=D1(l);
    S(m)=S1(l);
end
fid=fopen(fno,'w');
fprintf(fid,'频率(HZ)  阻尼(%%)  振型系数\n');
for k=1:mn
    fprintf(fid,'%10.4f  %10.4f  %10.6f\n',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);

土木年华

希望对大家有用，本人也在学习这一方法 希望与有共同方向的人共同商讨一下

ChaChing

建议先稍微交代内容及应用, 并在程序中加些注解, 这样好像较有完整性!
还有好像没给齐, 或许应该举个例子!?
最重要的原创或转贴应说明:@)