高斯白噪声的循环谱

N=2048;fs=1;df=fs/64;dalpha=fs/N;noise=wgn(1,2048,0);[Sx,alphao,fo]=fam(source_data2,fs,df,dalpha)ｆａｍ函数如下：function[Sx,alphao,fo]=fam(x,fs,df,dalpha)%************************************************%x为输入信号向量%fs为采样率%df为频率分辨率%dalpha为循环频率分辨率%注意：需满足dfdalpha才能得到较高的效果，fs/dalpha为进行采样计算功率谱的码元个数%%fam算法比ssca(autossca)算法快,但效果一样.都是求循环平稳与谱%************************************************%ifnargin4|nargin4%error('Wrongnumberofarguments.');%end%-----------DefinitionofParameters-----------%Np=pow2(nextpow2(fs/df));%输入信道数L=Np/4;%同一行连续列的相邻点的偏移P=pow2(nextpow2(fs/dalpha/L));%信道矩阵的列数N=P*L;%输入数据的点数%----------InputChannelization----------------%iflength(x)Nx(N)=0;elseiflength(x)Nx=x(1:N);end%m1=1:Fs:N;%x=(m1);%length(x);NN=(P-1)*L+Np;xx=x;xx(NN)=0;xx=xx(:);X=zeros(Np,P);fork=0:P-1X(:,k+1)=xx(k*L+1:k*L+Np);%输入数据xx的1到Np存入X第一列，L+1到L+Np存入第二列，2L+1到2L+Np存入第三列，以此类推。即L为偏移，Np为每段长度。end%------------------Windowing------------------%a=hamming(Np);XW=diag(a)*X;%每段加窗%XW=X;%不加窗，可与加窗的效果作比较%----------------FirstFFT--------------------%XF1=fft(XW);clearXW;%%%?????XF1=fftshift(XF1);XF1=[XF1(:,P/2+1:P)XF1(:,1:P/2)];%这两行的功能是把傅里叶变换的结果上半块和下半块交换,左半块和右半块互换%-------------Downconversion------------------%E=zeros(Np,P);fork=-Np/2:Np/2-1form=0:P-1E(k+Np/2+1,m+1)=exp(-i*2*pi*k*m*L/Np);endendXD=XF1.*E;clearXF1;%%%%???XD=conj(XD');%XD转置的复共轭clear('XF1','E','XW','X','x');%-----------Multiplication--------------------%XM=zeros(P,Np^2);fork=1:Npforq=1:NpXM(:,(k-1)*Np+q)=(XD(:,k).*conj(XD(:,q)));endendclearXD;%------------SecondFFT-----------------------%XF2=fft(XM);%clearXM;XF2=fftshift(XF2);XF2=[XF2(:,Np^2/2+1:Np^2)XF2(:,1:Np^2/2)];%length(XF2);XF2=XF2(P/4:3*P/4,:);M=abs(XF2);clearXF2;%Absolutevalueandcomplexmagnitude%%%%%%%%%%%%%%%%频率分辨率和循环频率分辨率alphao=(-1:1/N:1)*fs;%aboutthevariableN!!!fo=(-0.5:1/Np:0.5)*fs;Sx=zeros(Np+1,2*N+1);%aboutthevariableN!!!fork1=1:P/2+1fork2=1:Np^2ifrem(k2,Np)==0q=Np/2-1;elseq=rem(k2,Np)-Np/2-1;endk=ceil(k2/Np)-Np/2-1;p=k1-P/4-1;alpha=(k-q)/Np+(p-1)/L/P;f=(k+q)/2/Np;ifalpha-1|alpha1k2=k2+1;elseiff-0.5|f0.5k2=k2+1;elseifrem(k+q,2)==0&&rem(1+N*(alpha+1),1)==0kk=1+Np*(f+0.5);qq=1+N*(alpha+1);%aboutthevariableN!!!Sx(kk,qq)=M(k1,k2);endendenda=max(max(Sx));%归一化Sx=Sx./a;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下为算法画图figure;title('SCDestimateusingFAM')subplot(2,2,1)mesh(alphao,fo,Sx);xlabel('Cyclefrequency(alpha)')('frequency(f)')zlabel('Sx')gridon%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subplot(2,2,2)contour(alphao,fo,Sx)xlabel('Cyclefrequency(alpha)')ylabel('frequency(f)')gridon%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subplot(2,2,3)plot(alphao,10*log10((Sx(Np/2+1,:))))%f=0的情况未归一化,两边都有.xlabel('Cyclefrequency(alpha)')title('f=0')gridon%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subplot(2,2,4)plot(fo,10*log10(Sx(:,N+1)))%alpha=0的情况('frequency(f)')title('alpha=0')gridon;title('CMMBSignal');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(2);%title('SCDestimateusingFAM')%mesh(alphao,fo,Sx);%xlabel('Cyclefrequency(alpha)')%ylabel('frequency(f)')%zlabel('Sx')%title('CyclicSpectrumofCMMBsignal');%gridon;