prony.m

这是一个 York大学开发的FDTD代码。很好用的工具

function [C,Sp]=prony(f,Ts,Np)
% [C,Sp]=prony(f,Ts,Np)
% Use the method of Prony to determine the representation of a time response
% as a set of exponential basis functions such that:
%    f(N*Ts)=C0*exp(Sp0*N*Ts) + C1*exp(Sp1*N*Ts) + ... + CNp-1*exp(SpNp-1*N*Ts) 
% Need to check transpose v conjugate transpose in here !
%  J. F. Dawson 27 Oct 1999

Nt=size(f,1);

%Make F and fdot matrix of samples
F=f(1:Nt-Np);
for cnt=2:Np
    F=[F,f(cnt:Nt-Np+cnt-1)];
end

fdot=f(Np+1:Nt);

%size(F)
%size(fdot)

%Make and B matrices
A=F'*F;
B=-F'*fdot;

%
disp('A: Rank,Norm,Det')
rank(A)
norm(A)
det(A)
disp('-------------')


% Solve alpha=inv(A)*B Compute alpha coefficients 
alpha=A \ B;

%Check by back substitution
%[F*alpha+fdot] %should be zero


% Now compute the poles by finding the roots of the alpha polynomial
% mu^Np + alpha_Np-1*mu^Np-1 + alpha_Np-2*mu^Np-2 .... alpha_0
% =(mu - mu_0)(mu - mu_1)....(mu - mu_Np-1) = 0
% we want the mu_i values where mu_i=exp(Sp_i*Ts)

% find roots converting alpha array to the form required for root()
mu=roots([1;alpha(Np:-1:1)]);

%Frequency/time constant of each pole
Sp=log(mu)./Ts;


%Now we can compute the amplitudes coefficients:
mu=mu.';
MU=ones(1,Np);
MU=[MU;mu];
for cnt=2:Nt-1
    MU=[MU;mu.^cnt];
end

D=MU'*MU;
E=MU'*f;
C=D\E;

%Return results

%---------- Fix it for now.... --------------------

%Amplitude coefficients for each pole
%C=[1;0.5];

%Frequency/time constant of each pole
%Sp=[-0.01;-1];