# use prony method for damped exponentials

let us consider represented of damping exponential model by prony method,there is source code

y=zeros(1,N);
for i=1:N
y(i)=x(800*i);
end
d=zeros(1,N/2);
for i=1:N/2
d(i)=y(i+N/2);
end
D=zeros(N/2,N/2);
for i=1:N/2
for j=N/2:-1:1
D(i,-j+N/2+1)=y(i+j-1);
end
end
a=pinv(D)*d';

muhat=roots([1,-a']);
U=zeros(N,N/2);
for i=1:N
for j=1:N/2
U(i,j)=muhat(j,1)^(i-1);
end
end
C=pinv(U)*y';


and equation of model is following with there solving procedures : with their description : and solving strategies : Find Roots of charactreristic polynomial formed from the linear prediction coefficients Solve the original set of linear equations to yield the estimates of the exponential amplitude and sinusoidal phase i want following thing: we have given signal and i want input be it's length and L,or function should have following form

function [amplitudes,damping_factor,phase,frequency]=prony(y,n,L)

%n-length(y)
%L-number of complex exponential,


how can i continue?how to change given code in my case?thanks in advance

you have C and U, and C can be written as C_i = (A_i/2)*exp(j*phase) U_i can be written as exp(damp_i + j*2*pi*f_i)*T its should be easy to extract the info from that. T is sampling period.

%to get damping factor i and frequency i use something like this
%assume damping factor at sample 1 is 0.3, frequency at 1 is 20
%and out sampling frequency is 500 Hz
damp1 = [0.3 0.4 0.5];
f1 = [20 30 40];
T = 1/500;
%assume this is your ui vector from th ealgorithm
ui = exp((damp1 + j*2*pi*f1)*T);
%get the damping and frequencies
lv1 = log(ui);
r_damp1 = real(lv1)/T;
r_f1 =imag(lv1)/(T*2*pi);

%assume this is your Ci vector from algorithm
Ai = [0.1 0.2 0.3];
phasei=[pi/8 pi/4 pi/2];
Ci = (Ai/2).*exp(j*phasei);
%extract phase and amplitude
r_Amp = abs(Ci)*2;
r_phase = angle(Ci);

• how can i know T by the way? Apr 7 '14 at 12:09
• whatever your input signal to the algorithm was sampled at. these things are usually known before any algorithm processing. for audio its usually 44.1KHz or 48KHz sampling frequency, for speech its lower. Apr 7 '14 at 12:14
• i have outputs,C and U,do i need T now?or simple i can continue? Apr 7 '14 at 12:14
• what is your input signal ? is it speech, audio ? Apr 7 '14 at 12:30
• damp = real(log(U))/T, freqs=imag(log(U))/(T*2*pi),Amp = abs(C)*2,phase = angle(C) Apr 7 '14 at 19:52