pisar.m

来自「谱估计 (1)修正协方差法 (2)多重信号分类（MUSIC）算法 (3)E」· M 代码 · 共 72 行

function [A,F]=pisar(xt,sin_num)

%   PISAR Pisarenko Harmonic Decomposition
%   
%   PISAR(XT,SIN_NUM) Subplots the input signal (Xt) and the
%   signal generated by the Pisarenko Harmonic Decomposition 
%   algorithm for a sum of SIN_NUM sinusoid waves.
%   
%   [A F]=PISAR(XT,SIN_NUM) Generates SIN_NUM sinusoids, where the A
%   column-vector represents the amplitudes and the column-vector F
%   represents the normalized frequencies of those sinusoids.
%   
%   Considerations:
%      XT must be a column vector;
%      For best results --  SIN_NUM << length of XT
%
%   Plese send any comments to [e-mail]
%   			       Luis E. Palafox Maestre
%

[N t]=size(xt);
rxx=xcorr(xt,'biased');
rxx=rxx(N:(2*sin_num)+N);

%Frequencies estimation
Rxx=toeplitz(rxx);
ev=eig(Rxx);
[S i]=min(ev);
[V D]=eig(Rxx);
a=V(:,i);
rts=roots(a);
w_est=[];
for i=1:sin_num
   w_est(i)= abs(angle(rts(2*i)));
end
F=(w_est/(2*pi))';

%Amplitudes estimation
mcos=[];
for n=1:sin_num
   vcos=[];
   for i=1:sin_num
      vcos=[vcos cos(n*w_est(i))];
   end
   mcos=[mcos; vcos];
end
rxx=rxx(2:sin_num+1);
rxx=2*rxx;
A=inv(mcos)*rxx;
A=A^(1/2);

if nargout==0,
  xe=[];
  for n=1:N
    xe(n)=0;
    for i=1:sin_num
      xe(n)= xe(n)+(A(i)*cos(w_est(i)*n));
    end
  end
  f=figure;
  subplot(2,1,1);
  plot(xt);
  title('Input Signal')
  xlabel('n');
  ylabel('x(n)');
  subplot(2,1,2);
  plot(xe);
  title('Estimated Signal');
  xlabel('n');
  ylabel('x(n)');
end