高斯白噪聲的循環(huán)譜

1、精選優(yōu)質(zhì)文檔-傾情為你奉上N=2048;fs=1;df=fs/64;dalpha=fs/N;noise=wgn(1,2048,0);Sx,alphao,fo = fam(source_data2,fs,df,dalpha)函數(shù)如下：function Sx,alphao,fo = fam(x,fs,df,dalpha)%*% x為輸入信號(hào)向量% fs為采樣率% df為頻率分辨率% dalpha為循環(huán)頻率分辨率% 注意：需滿足df>>dalpha才能得到較高的效果，fs/dalpha為進(jìn)行采樣計(jì)算功率譜的碼元個(gè)數(shù)%fam算法比ssca(autossca)算法快,但效果一樣.都是求循環(huán)平

2、穩(wěn)與譜%*% if nargin > 4 | nargin < 4%     error('Wrong number of arguments.');% end%-Definition of Parameters-%Np = pow2(nextpow2(fs/df);     %輸入信道數(shù)L = Np/4;                     &

3、#160; %同一行連續(xù)列的相鄰點(diǎn)的偏移P = pow2(nextpow2(fs/dalpha/L);%信道矩陣的列數(shù)N = P*L;                        %輸入數(shù)據(jù)的點(diǎn)數(shù)%-Input Channelization-%if length(x) < N    x(N) = 0;elseif length(x) > N    x = x(1:N);en

4、d% m1=1:Fs:N;% x=(m1);% length(x);NN = (P-1)*L+Np;xx = x;xx(NN) = 0;xx = xx(:);X = zeros(Np,P);for k = 0:P-1    X(:,k+1) = xx(k*L+1:k*L+Np);      %輸入數(shù)據(jù)xx的1到Np存入X第一列，L+1到L+Np存入第二列，2L+1到2L+Np存入第三列，以此類推。即L為偏移，Np為每段長(zhǎng)度。end%-Windowing-%a = hamming(Np);XW = diag(a)*X;  &#

5、160;%每段加窗% XW = X;         %不加窗，可與加窗的效果作比較%-First FFT-%XF1 = fft(XW);  clear XW;    %  ?XF1 = fftshift(XF1);   XF1 = XF1(:,P/2+1:P) XF1(:,1:P/2);  %這兩行的功能是把傅里葉變換的結(jié)果上半塊和下半塊交換,左半塊和右半塊互換%-Downconversion-%E = zeros(Np,P);for

6、 k = -Np/2:Np/2-1    for m = 0:P-1        E(k+Np/2+1,m+1) = exp(-i*2*pi*k*m*L/Np);    endendXD = XF1.*E;   clear XF1;  %  ?XD = conj(XD');                  

7、    %XD轉(zhuǎn)置的復(fù)共軛clear ('XF1', 'E', 'XW', 'X', 'x'); %-Multiplication-%XM = zeros(P,Np2);for k = 1:Np    for q = 1:Np        XM(:,(k-1)*Np+q) = (XD(:,k).*conj(XD(:,q);    endendclear XD;%-Second FFT-

8、%XF2 = fft(XM);  % clear XM;XF2 = fftshift(XF2);XF2 = XF2(:,Np2/2+1:Np2) XF2(:,1:Np2/2);% length(XF2);XF2 = XF2(P/4:3*P/4,:);M = abs(XF2);   clear XF2;  %Absolute value and complex magnitude%頻率分辨率和循環(huán)頻率分辨率alphao = (-1:1/N:1)*fs;     %about the variabl

9、e N!fo = (-0.5:1/Np:0.5)*fs;Sx = zeros(Np+1,2*N+1);     %about the variable N!for k1 = 1:P/2+1    for k2 = 1:Np2        if rem(k2,Np) = 0            q = Np/2-1;         &

10、#160;  else            q = rem(k2,Np)-Np/2-1;         end        k = ceil(k2/Np)-Np/2-1;         p = k1-P/4-1;        al

11、pha = (k-q)/Np+(p-1)/L/P;        f = (k+q)/2/Np;        if alpha < -1 | alpha > 1            k2 = k2+1;        elseif f < -0.5 | f > 0.5     &#

12、160;      k2 = k2+1;        elseif rem(k+q,2)=0 && rem(1+N*(alpha+1),1)=0            kk = 1+Np*(f+0.5);            qq = 1+N*(alpha+1);   %about the vari

13、able N!            Sx(kk,qq) = M(k1,k2);                end    endenda=max(max(Sx);     %歸一化Sx=Sx./a;%以下為算法畫圖figure;title('SCD estimate using FAM')subplot(2,2,1)m

14、esh(alphao,fo,Sx);xlabel('Cycle frequency(alpha)')ylabel('frequency(f)')zlabel('Sx')grid on%subplot(2,2,2)contour(alphao,fo,Sx)xlabel('Cycle frequency(alpha)')ylabel('frequency(f)')grid on%subplot(2,2,3)plot(alphao,10*log10(Sx(Np/2+1,:)%f=0的情況未歸一化,兩邊都有.xlabel('Cycle frequency(alpha)')title('f=0')grid on% subplot(2,2,4)plot(fo,10*log10(Sx(:,N+1)%alpha=0的情況xlabel('frequency(f)')title('alpha=0')grid on;title('CMMB Signal');% figur