lyj.m

一种基于四阶累积量的谐波恢复的方法

A1=0;A2=0;F1=0;F2=0;
for q=1:10
r=unifrnd(-pi,pi,1,2);%
a1=0.6;a2=1;f1=800;f2=2400;fs=8000;s1=0;s2=0;
n=1:102400;
f1=f1/fs;f2=f2/fs;
x(n)=a1*exp(j*(2*pi*f1*n+r(1)))+a2*exp(j*(2*pi*f2*n+r(2)));
h=[0.021496113113274748 0.039135016469765589 0.003128535040244775 -0.11993091620746835 0.2666238654318458 0.66998625441928483 0.2666238654318458 -0.11993091620746835 0.003128535040244775 0.039135016469765589 -0.021496113113274748   ];
%lowpass fpass:3000 fstop:5000
%n=wgn(1,102390,-40)+j*wgn(1,102390,-40);
%n=unifrnd(0,0.0125,1,1014)+j*unifrnd(0,0.0125,1,1014);
n=random('poisson',0.000125,[1,102390])+j*random('poisson',0.000125,[1,102390]);%瑞利分
%布
n=conv(n,h);
for i=1:102400
    s1=s1+x(i)*conj(x(i));
    s2=s2+n(i)*conj(n(i));
end
SNR=10*log10(abs(s1)/abs(s2));%计算信噪比


y=x+n;
%M=4;
%cum_4y=cumest(y,4,M,length(y),0,'biased',0,0);
%cum_4y=cum44(y,M);
%cum_4y=(cum_4y(M+1:2*M+1)+conj(cum_4y(M+1:-1:1)))/2;
%y_cum4=cum_4y/cum_4y(1);


%v=[y_cum4(3),y_cum4(2),y_cum4(1);y_cum4(4),y_cum4(3),y_cum4(2);y_cum4(5),y_cum4(4),y_cum4(3)];
%v=inv(v);
%a=[0.000000000000000000001;0.000000000000000000000001;0.000000000000000000000000000001];
%a=v*a;
M=2;
cum_4y=cumest(y,4,M,length(y),0,'biased',0,0);
%cum_4y=cum44(y,M);
cum_4y=(cum_4y(M+1:2*M+1)+conj(cum_4y(M+1:-1:1)))/2;
y_cum4=cum_4y/cum_4y(1);
Rxx=toeplitz(y_cum4,conj(y_cum4));

ev=eig(Rxx);
[S i]=min(ev);
[V D]=eig(Rxx);
a=V(:,i);
rts=roots(a);
w_est=[];
for i=1:2
   w_est(i)= angle(rts(i));%估计频率
end
F=(w_est/(2*pi))'*8000;
if(F(1)>1500)
    x=F(1);
    F(1)=F(2)
    F(2)=x;
end

F1=F(1)+F1;
F2=F(2)+F2;
%figure
  %stem(F,ones(1,length(F)));
  %title('Frequencies estimated by cum4 method');
m=[1,1; exp(j*2*pi*F(1)/fs),exp(j*2*pi*F(2)/fs)];%估计幅度
G=inv(m);
h=y_cum4(2:3);
t= (-1)*G*h;
for i=1:2
Amat(i)=sqrt(sqrt(abs(t(i))));
end
if(Amat(1)>0.8)
     x=Amat(1);
    Amat(1)=Amat(2)
    Amat(2)=x;
end
A1=A1+Amat(1);
A2=A2+Amat(2);


end
F1=F1/10;
F2=F2/10;
A1=A1/10;
A2=A2/10;