MATLAB谐波信号双谱估计

clcclearcloseallN=128*64;n=[0:N-1];Nfft=128;Fs=1;t=n/Fs;f1=0.6381/2/pi;f2=0.8345/2/pi;f3=f1+f2;f4=0.4909/2/pi;f5=1.7671/2/pi;f6=f4+f5;fai1=pi/6;fai2=pi/3;fai3=fai1+fai2;fai4=5*pi/4;fai5=6*pi/4;fai6=pi;s1=cos(2*pi*f1*t+fai1);s2=cos(2*pi*f2*t+fai2);s3=cos(2*pi*f3*t+fai3);s4=cos(2*pi*f4*t+fai4);s5=cos(2*pi*f5*t+fai5);s6=cos(2*pi*f6*t+fai6);s=s1+s2+s3+s4+s5+s6;%s=awgn(s,0,'measured');noi=zeros(1,N);%高斯白噪声noi=randn(1,N);s=s+noi;%有色噪声nys=0;if(nys==1)Ln=length(s);Nd=[1-0.50.70.1];Nc=[10.50.2];nd=length(Nd)-1;nc=length(Nc)-1;xik=zeros(nc,1);ek=zeros(nd,1);xi=randn(Ln,1);forkn=1:Lne(kn)=-Nd(2:nd+1)*ek+Nc*[xi(kn);xik];fornn=nd:-1:2ek(nn)=ek(nn-1);endek(1)=e(kn);fornn=nc:-1:2xik(nn)=xik(nn-1);endxik(1)=xi(kn);ends=s+0.25*e;noi=0.25*e;end%信噪比S_N=10*log10(sum(abs(s).^2)/sum(abs(noi).^2));[bspec,waxis]=bispecd(s,Nfft,1,64,32);closeall%双谱信号（包含相干信号的幅度和相位信息，幅度受信号的功率谱影响）SS=Nfft/2+1;waxis1=waxis(SS:end)*Fs;figuresubplot(121)bspec1=bspec(SS:end,SS:end);contour(waxis1,waxis1,abs(bspec1),4);gridontitle('双谱估计二维等高线图')xlabel('f1/Hz'),ylabel('f2/Hz')[X1,Y1]=meshgrid(waxis1,waxis1);subplot(122)mesh(X1,Y1,abs(bspec1))axistightxlabel('f1/Hz')ylabel('f2/Hz')zlabel('幅度')title('双谱估计三维图');figureplot(t,s);axis([0,200,-8,8])xlabel('t')ylabel('幅度')title('时域信号')P=10*log10(abs(fft(s-mean(s),Nfft).^2)/(Nfft+1));f=(0:length(P)-1)*Fs/(length(P));figure()plot(f(1:Nfft/2),P(1:Nfft/2));title('功率谱估计');xlabel('f/Hz')ylabel('P/dB')%相干双谱信号（频率f1,f2非线性相位耦合产生的能量在f1+f2处总能量中所占的比例）mask=hankel([1:Nfft],[Nfft,1:Nfft-1]);%thehankelmask(faster)P_d=abs(fft(s,Nfft).^2)/(Nfft+1);P_d=P_d';bspec_d=zeros(Nfft,Nfft);Bspec_d=zeros(Nfft,Nfft);bspec_d=(P_d*P_d').*reshape(P_d(mask),Nfft,Nfft);Bspec_d=abs(bspec)./sqrt(bspec_d);figure()subplot(121)contour(waxis1,waxis1,abs(Bspec_d(SS:end,SS:end)),4);gridontitle('双相干谱估计二维等高线图')xlabel('f1/Hz'),ylabel('f2/Hz')subplot(122)mesh(X1,Y1,abs(Bspec_d(SS:end,SS:end)));axistightxlabel('f1/Hz')ylabel('f2/Hz')zlabel('幅度')title('双相干谱估计三维图');figure()subplot(211)diagBspec=diag(fliplr(bspec));plot(waxis1,abs(diagBspec(SS:end)));gridon;xlabel('f/Hz')ylabel('幅度')title('双谱估计切片图')subplot(212)diagBspec_d=diag(fliplr(Bspec_d));plot(waxis1,abs(diagBspec_d(SS:end)));gridon;xlabel('f/Hz')ylabel('幅度')title('双相干谱估计切片图')