feat_extr.rar_.ana_extr_extraction_feature extraction_feature_ex

function [gamma_max,sigma_ap,sigma_dp,P]=feat_extr(x,fc,fs,fx,at) %usage:extract key feature % gamma_max:apmplitude information % sigma_ap:absolute phase information % sigma_dp:direct phase information % P:spectrum symmetry %input parameters % x:signal received % fc:carry frequency % fs:sampple frequency % fx:message signal bandwidth,ie 8 khz % at:threshold for {a(i)} below which the estimation of the instantaneous % phase is sensitive to noise Ns=length(x); Xc=fft(x); fcn=fc*Ns/fs; %计算频谱对称性参数 p Pl=0;Pu=0; % for i=1:fcn % Pl=Pl+Xc(i)*Xc(i)'; % Pu=Pu+Xc(i+fcn+1)*Xc(i+fcn+1)'; % end; Pl=Xc(1:fcn)'*Xc(1:fcn); Pu=Xc((fcn+2):(2*fcn+1))'*Xc((fcn+2):(2*fcn+1)); P=(Pl-Pu)/(Pl+Pu); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % i=1:Ns; %正交下变频求hilbert变换 % I=cos(2*pi*fc*i/fs); % I=I'; % Q=sin(2*pi*fc*i/fs); % Q=Q'; % coef=fir1(30,fx/(fs/2)); % xI=x.*I; % xQ=x.*Q; % % s=fftfilt(coef,xI)+j*fftfilt(coef,xQ); %hilbert transform of x % s=xI+j*xQ; % a=abs(s); %instantaneous amplitude % fai=angle(s); %instantaneous phase % c=zeros(Ns,1); % if fai(2)-fai(1)>pi % c(1)=-2*pi; % else if fai(1)-fai(2)>pi % c(1)=2*pi; % end % end % for t=2:Ns-1 % if fai(t+1)-fai(t)>pi % c(t)=c(t-1)-2*pi; % else if fai(t)-fai(t+1)>pi % c(t)=c(t-1)+2*pi; % end % c(t)=c(t-1); % end % end % c(Ns)=c(Ns-1); % % fai_uw=fai+c(t); % t=1:Ns; % temp=2*pi*fc*t/fs; % temp=temp.'; % fai_nl=fai_uw-temp; %nonlinear instantaneous phase % fai_nl=fai_nl-mean(fai_nl); %centered % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z=hilbert(x); %计算gamma_max a=abs(z); ma=mean(a); an=a/ma; acn=an-1; Acn=fft(acn); gamma_max=( max(abs(Acn)) )^2/Ns; % y=imag(z); % 计算sigma_ap,sigma_dp % fai=zeros(Ns,1); % for t=1:Ns % if x(t)>0&y(t)>=0 fai(t)=atan(y(t)/x(t)); % else if x(t)<0&y(t)>=0 fai(t)=pi-atan(y(t)/x(t)); % else if x(t)<0&y(t)<=0 fai(t)=pi+atan(y(t)/x(t)); % else if x(t)>0&y(t)<=0 fai(t)=2*pi-atan(y(t)/x(t)); % else if x(t)==0&y(t)>0 fai(t)=0.5*pi; % else if x(t)==0&y(t)<0 fai(t)=1.5*pi; % end % end % end % end % end % end % end fai=angle(z); % fai=phase(z); fai_uw=unwrap(fai); t=1:Ns; temp=2*pi*fc*t/fs; fai_nl=fai_uw-(temp)'; fai_nl=fai_nl-mean(fai_nl);%centred t1=0;t2=0;t3=0;C=0; ipos=find(an>at); C=length(ipos); t1=sum(fai_nl(ipos).^2); t2=sum(abs(fai_nl(ipos))); t3=sum(fai_nl(ipos)); % t11=0;t21=0;t31=0; % for i=1:Ns % if(an(i)>at) % C=C+1; % t11=t11+fai_nl(i)^2; % t21=t21+abs(fai_nl(i)); % t31=t31+fai_nl(i); % end; % end; sigma_ap=sqrt(t1/C-(t2/C)^2); sigma_dp=sqrt(t1/C-(t3/C)^2);