function [poles,coeff,ampl,radii,freq,phase]=prony(x,N) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function [poles,coeff,ampl,radii,freq,phase]=prony(x,N) % % PRONY MODELING % % OUTPUT % f = frequencies of poles % r = radii of poles % c = coefficients % p = phases % % INPUT % x = original signal % N = number of poles to try % % Example: % n=0:15 % x=0.5*sin(2*pi*n/12)+(0.99.^n).*cos(2*pi*n/6.2)+0.97.^n % [poles,coeff,ampl,radii,freq,phase]=prony(x,5) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % B. Hutchins 1994 2001 2006 Cornell Univ. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clf % Compute Length 2N Signal Vector for Model % Extract x1 from full x idx=1:2*N; x1=x(1:2*N); % Look at input % Plot Input x and selected samples x1 as * figure(1) plot(x) hold on plot(idx,x1,'*') xrng=(max(x)-min(x))/10; axis([0 length(x)+1 min(x)-xrng max(x)+xrng]); hold off title('Input, and Samples to be Used *') % PRONY METHOD PART 1, Find the Poles % giving frequency and decay constants % % Solve for Difference Equation Coefficients % Set Up Matrix b and Vector d1 for K=1:N d1(K)=x1(N+K); for L = 1:N b(K,L)=x1(K+L-1); end end % Solve for Coefficients a a = inv(b)*d1'; % Find Poles d=[-a',1]; d=fliplr(d); % denominator of network = char. poly. poles=roots(d); % poles r=abs(poles); % pole radii an=angle(poles); % pole angles figure(2) % Plot poles p1=roots(d); pmax=max([abs(real(p1)); abs(imag(p1))]); sc=ceil(pmax)+0.1; n=0:500; r1=exp(j*2*pi*n/500); plot(r1,'g') grid hold on plot(real(p1),imag(p1),'x') hold off axis('equal') axis([-sc sc -sc sc]) % end plotting of poles % f=(an/(2*pi)); % actual pole frequencies % END PART 1 % % % % PRONY METHOD PART 2, Apply Initial Conditions % to find amplitude and phase for K=1:N for L = 1:N cc(K,L)=poles(L)^(K-1); end end for K=1:N g(K)=x1(K); end c=inv(cc)*g'; magc=abs(c); % END PART 2 % Test Reconstruction to length 10N for K=1:10*N xr(K)=0; for L=1:N xr(K)=xr(K)+c(L)*poles(L)^(K-1); end end xr=real(xr); % PLOTS figure(3) plot(x,'g') hold on pv1=[1:2*N]; xr1=xr(1:2*N); plot(pv1,xr1,'o') pv2=[(1+2*N):(1+10*N-1)]; xr2=xr((2*N+1):length(xr)); plot(pv2,xr2,'+') lpv1=length(pv1); pv2p=[pv1(lpv1), pv2]; xr2p=[xr1(lpv1), xr2]; plot(pv2p,xr2p,'r'); plot([0 10*N],[0 0],'c') mx=max([xr2 x]); mn=min([xr2 x]); rg=0.1*(mx-mn); axis([0 length(xr2)+5 mn-rg mx+rg]); hold off % phase relative to a cosine input p=360*angle(c)/(2*pi); radii=r; coeff=c; ampl=2*abs(coeff); % correct amplitude for real poles for k=1:N if abs(imag(poles(k)))<0.000001 ampl(k)=ampl(k)/2; end end freq=f; phase=p;