prony算法matlab编程?

x=0:0.05:10% 原理想信号

y=220*exp(-0.25*x).*cos(2*pi*1.5*x+pi)+110*exp(0.25*x).*cos(2*pi*1.0*x+pi/2)

p=6% p 为模型阶数

sf=1/0.05 % sf 为采样频率

[F,D,A,theta]=prony_my(y,p,sf)

% 调用prony 算法函数

function [F,D,A,theta]=prony_my(y,p,sf)

nm=2 * p

n=fix(length(y)/2)

dt=1/sf% 由采样频率宴辩态计算时间间隔

L=length(h)

M=L/2

for k= 1: nm

x1(:,k) =h(k: M-1+k)

end

for k=1: M

x2(k,:) =-h(nm+k)

end

B=x1\x2% 最小二乘法求解prony 多项式系数

B(nm+1) =1

B1=B(nm+1 : -1: 1)

V=roots(B1)% 求出特征方程的根

F1=abs(log(V)) /晌源 (2*pi*dt)%初步计算频率

D1=log(abs(V))/dt% 初步计算衰减因子

for k=0: (2*n-1)

Va(k+1,:) = conj(V').^k%生成范得蒙矩阵Va

end

S1=2*(inv(conj(Va') * Va) *conj(Va') * h)

[F2, I] =sort(F1)% 将频率F1从小到大排序

m=0

for k=1: nm-1

if F2(k) ~=F2(k+1) % 消除其中的共灶搭轭项

continue

end

m=m+1

ii=I(k)% 提取有效的元素序数

F(m) =F1(ii)

D(m) =D1(ii)

A(m) =abs(S1(ii))

theta(m)=angle(S1(ii))

end

function [b,a] = prony(h, nb ,na)

%PRONY Prony's method for time-domain IIR filter design.

% [B,A] = PRONY(H, NB, NA) finds a filter with numerator order

% NB, denominator order NA, and having the impulse response in

% vector H. The IIR filter coefficients are returned in

% length NB+1 and NA+1 row vectors B and A, ordered in

% descending powers of Z. H may be real or complex.

%

% If the largest order specified is greater than the length of H,

% H is padded with zeros.

%

% See also STMCB, LPC, BUTTER, CHEBY1, CHEBY2, ELLIP, INVFREQZ.

% Author(s): L. Shure, 5-17-88

% L. Shure, 12-17-90, revised

% Copyright 1988-2002 The MathWorks, Inc.

% $Revision: 1.7 $ $Date: 2002/03/28 17:30:10 $

% References:

% [1] T.W. Parks and C.S. Burrus, Digital Filter Design,

% John Wiley and Sons, 1987, p226.

K = length(h) - 1

M = nbN = na

if K <= max(M,N) % zero-pad input if necessary

K = max(M,N)+1

h(K+1) = 0

end

c = h(1)

if c==0% avoid divide by zero

c=1

end

H = toeplitz(h/c,[1 zeros(1,K)])

% K+1 by N+1

if (K >N)

H(:,(N+2):(K+1)) = []

end

% Partition H matrix

H1 = H(1:(M+1),:)% M+1 by N+1

h1 = H((M+2):(K+1),1)% K-M by 1

H2 = H((M+2):(K+1),2:(N+1))% K-M by N

a = [1-H2\h1].'

b = c*a*H1.'

prony_my函数部分应该单独保存为prony_my.m文件，也就是从function开头的那一句开兄绝始应该专门保存毕尘段为prony_my.m文件。然后在命令窗口输入：

x=0:0.05:10% 原理想信号

y=220*exp(-0.25*x).*cos(2*pi*1.5*x+pi)+110*exp(0.25*x).*cos(2*pi*1.0*x+pi/2)

p=6% p 为模型阶数

sf=1/0.05 % sf 为采样频率

[F,D,A,theta]=prony_my(y,p,sf)

% 调用prony 算法函手誉数

---