www.gusucode.com > 生物信号工具箱 > 生物信号工具箱/生物信号工具箱/Biosignal/Biosignal/Adaptive_Algorithms/harmonic_recon.m

    function [varargout]=harmonic_recon(varargin)
%
%[xhat,mse,theta,theta_curve,b,a,amp,err]=harmonic_recon(x,Fs,P,theta,mu,r,M)
%
%
% Reconstructs a signal based on it's harmonic components estimated through the
% harmonic_est function. Arguments are:
%
% Input:
%
% x        - (Nx1 required) signal to be estimated.
% Fs       - (1x1 required) Sampling frequency (Hz).
% P        - (1x1 required) Number of initial harmonics to used to estimate the
%            fundamental frequency. Once the fundamental is estimated the
%            signal is reconstructed with harmonics up to the nyquist rate.
% theta    - (1x1 optional) Initial estimate of the fundamental frequency (Hz).
% mu       - (1x1 optional) Step size for the harmonic_est function.
% r        - (1x1 optional) Controls the 'width' of the nulls used in the harmonic_est
%             function.
% M        - (1x1 optional) Scalar number of samples from which to compute xhat
%            (default is N, same size as x).
% Outputs:
%
% xhat          - (Nx1 required) Reconstructed signal from harmonics.
% mse           - (1x1 optional) Mean square error.
% theta         - (1x1 optional) Final fundamental frequency estimate.
% theta_curve   - (Nx1 optional) Learning curve of the fundamental estimate (useful
%                 for checking convergence of the algorithm).
% b             - (Pmx1 optional) Comb filter MA coefficients.
% a             - (Pmx1 optional) Comb filter AR coefficients.
% amp           - (Pmx1 optional) Amplitude of harmonic components
%
%
%
% % Example
% load ecg_ex
% theta=1.8*pi*2/Fs;
% tm=[0:1/Fs:(length(x)-1)/Fs]';
% mu=10^-10;
% r=0.9;
% P=25;
% N=length(x);
% M=1.5*N;
% [xhat,mse,theta,theta_curve,b,a,amp]=harmonic_recon(x,Fs,P,theta,mu,r,M);
%
%
% figure
% subplot(211)
% plot(xhat)
% hold on
% plot(x,'r')
% legend('Recon','True Sig')
%
% subplot(212)
% title('Spectral Analysis')
% Ff=[1:length(amp)].*theta;
% stem(Ff,amp./max(amp),'b-o')
% hold on;grid on
% [Pxx,Ff]=pwelch(x,ones(length(x),1),0,[],Fs);
% plot(Ff,sqrt(Pxx)./max(sqrt(Pxx)),'r')
% legend('Estimated Harmonics', 'Estimated Signal Spectrum')
%
%
%
%
%  Writen by Ikaro Silva 2010 version 1.0.1
%

x=varargin{1};
[N,L]=size(x);
Fs=varargin{2};
P=varargin{3};
theta=[];
mu=10^-10;
r=0.9;
M=N;
if(nargin>=4 && ~isempty(varargin{4}))
    theta=varargin{4};
end
if(nargin>=5 && ~isempty(varargin{5}))
    mu=varargin{5};
end
if(nargin>=6 && ~isempty(varargin{6}))
    r=varargin{6};
end
if(nargin>=7 && ~isempty(varargin{7}))
    M=varargin{7};
end


%Estimate Harmonics
[theta,theta_curve,b,a,err]=harmonic_est(x,P,Fs,theta,mu,r);

mx=mean(x);
x=x-mx;
tm=[0:1/Fs:(M-1)/Fs]';
xhat=0;
Pm=floor(Fs/(2*theta));
y=zeros(M,1);
amp=zeros(Pm,1);
for p=1:Pm
    sig=cos(2*pi*tm(1:N)*p*theta);
    [R,lag]=xcorr(sig,x);
    [trash,dly]=max(R);
    phi=lag(dly)/Fs;
    sig=cos(2*pi*(tm+phi)*p*theta);
    amp(p)=sig(1:N)'*x./(sig(1:N)'*sig(1:N));
    sig=amp(p).*sig;
    y=y + sig;
end

vx=std(x);
vy=std(y);
y=y.*vx/vy;
y=y+mx;
varargout(1)={y};


if(nargout>1)
    mse=mean((y(1:N)-x(1:N)).^2);
    varargout(2)={mse};
end
if(nargout>2)
    varargout(3)={theta};
end
if(nargout>3)
    varargout(4)={theta_curve};
end
if(nargout>4)
    varargout(5)={b};
end
if(nargout>5)
    varargout(6)={a};
end
if(nargout>6)
    varargout(7)={amp};
end
if(nargout>7)
    varargout(8)={err};
end