dsp谱估计.m

分别使用自相关函数法和周期图的方法对信号进行功率谱估计

N = input('输入数据总长度N：');
M = input('输入自相关函数个数（矩形窗宽度）M：');
K = input('输入数据分组的段数K：');
Ns = input('输入正弦个数Ns：');
%各个正玄信号的数字频率
fn = zeros(1,Ns);
%各个正玄信号的相位
phi = zeros(1,Ns);
%各个正玄信号的幅度
A = zeros(1,Ns);
for i = 1:Ns,
    sprintf('第%d个正弦信号的幅度:',i)
    A(i) = input('');
    sprintf('第%d个正弦信号的数字频率:',i)
    fn(i) = input('');
    sprintf('第%d个正弦信号的相位:',i)
    phi(i) = input('');
end
sigmaW = input('输入复高斯白噪声信号的方差：');
%生成N个方差等于sigmaW的复高斯白噪声W(n)
W = wgn(1, N, sigmaW, 1, 'linear', 'complex');
%信号x(n)
x = zeros(1,N);
for i = 1:N,
    for j = 1:Ns,
        x(i) = x(i) + complex(A(j)*cos(fn(j)*i+phi(j)),(-1)*A(j)*sin(fn(j)*i+phi(j)));
    end
    x(i) = x(i) + W(i);
end
%信号的自相关函数Rxx--------它是偶函数
Rxx = zeros(1,2*M-1);
for m = 1:M,
    te=m+M-1;
    for i = 1:N-m,
        Rxx(te) = Rxx(te) + x(i)*conj(x(i+m-1));
    end    
    Rxx(te) = Rxx(te)/(N-m+1);
    Rxx(M-m+1)=Rxx(te);
end
%0~pi 之间的角频率
fSample = zeros(1,128);
for i = 1:128,
    fSample(i) = (i-1)/128*pi;
end
MRect = zeros(1,128);
MHamming = zeros(1,128);
for i = 1:128,
    for m = -M+1:M-1,
          te=m+M;
           MRect(i) = MRect(i) + Rxx(te)*complex(cos(-m*fSample(i)),sin(-m*fSample(i))); 
           MHamming(i) = MHamming(i) + (0.54 + 0.46*cos(pi*m/M))*Rxx(te)*complex(cos(-m*fSample(i)),sin(-m*fSample(i)));
           
    end
end
i = 1:128;
%plot(i,abs(MRect(i)),i,abs(MHamming(i)));  
plot(i,MRect(i));
title('矩形窗时功率谱');
figure;
plot(i,MRect(i),i,MHamming(i));
title('矩形窗和Hamming窗时的功率谱');
legend('矩形窗','Hamming窗',0);
%Bartlett平均周期图方法
NF = input('输入FFT点数NF：');
%由已知可得到每段长度L为：
L = N/K;
Xi = zeros(1,NF);
Xi_FFT = zeros(1,NF);
I_w = zeros(K,NF);
P_w = zeros(1,NF);
for i = 1:K,
    %每段中的x(n)----Xi(n);i--------第几段
    for j = 1:L,
        Xi(j) = x(j+(i-1)*L);
    end
    %算出L点的FFT,后又算出I(w)
    Xi_FFT = fft(Xi,NF);
    for j = 1:NF,
        %j表示频率分点即Wi=2*pi*(j-1)/NF,其中1<=j<=NF
        I_w(i,j) = (abs(Xi_FFT(j)))^2/L;        
    end   
end
for j = 1:NF,
    for i = 1:K,
        P_w(j) = P_w(j) + I_w(i,j)/K;
    end
end
figure;
plot(P_w);
title('周期图法－－矩形窗');
%改进Bartlett平滑平均周期图方法
%由已知可得到每段长度L为：
Xi = zeros(1,NF);
Xi_FFT = zeros(1,NF);
I_w = zeros(K,NF);
P_w = zeros(1,NF);
U = 0;
for j = 1:L,
    U = U + (0.54+0.46*cos(pi*(j-1)/L))^2;
end
U = U/L;
for i = 1:K,
    %每段中的x(n)----Xi(n);i--------第几段
    for j = 1:L,
        Xi(j) = x(j+(i-1)*L)*(0.54+0.46*cos(pi*(j-1)/L));
    end
    %算出L点的FFT,后又算出I(w)
    Xi_FFT = fft(Xi,NF);
    for j = 1:NF,
        %j表示频率分点即Wi=2*pi*(j-1)/NF,其中1<=j<=NF
        I_w(i,j) = (abs(Xi_FFT(j)))^2/(L*U);        
    end   
end
for j = 1:NF,
    for i = 1:K,
        P_w(j) = P_w(j) + I_w(i,j)/K;
    end
end
figure;
plot(P_w);
title('周期图法－－Hamming窗');