feat_extr.m

ana modulation feature extraction

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);