间接法求功率谱密度之自相关函数---基于MATLAB

方法：先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅里叶变换，便得到x(n)的功率谱估计。

（1）、自相关函数
求解步骤：（1）移位；（2）相乘；（3）相加。

function y = myXcorr(x)
%输入x：长度为n的一维行向量
%输出y：长度为2n-1的一维行向量
%求解步骤：(1)移位;(2)相乘;(3)相加.
n = length(x); %输入行向量的长度
y = zeros(1,2*n-1);  %提前分配输出行向量的空间
for i = 1:n
    %利用前后补0进行移位操作
    x1 = [zeros(1,n-i), x];  
    x2 = [x, zeros(1,n-i)];  
    y([i,2*n-i]) = x1*x2';   %进行相乘和相加操作
end

验证：

clc
clear

a1 = [8 7 2 4 6 8];
a2 = [5 6 8 7 4 1 2];

b1 = xcorr(a1)  %MATLAB内置函数计算结果
c1 = myXcorr(a1)  %自编MATLAB函数计算结果

b2 = xcorr(a2)  %MATLAB内置函数计算结果
c2 = myXcorr(a2)  %自编MATLAB函数计算结果

结果：

（2）求功率谱密度

clc;
clear;

Fs = 1000;  %采样频率
n = 0:1/Fs:1;

xn = cos(2*pi*40*n) + 3*cos(2*pi*100*n) + randn(size(n));   %含噪声信号
nfft = 1024;

% 间接法
%%% 间接法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅里叶变换，便得到x(n)的功率谱估计。

%计算序列的自相关函数
cxn = xcorr(xn);  
mycxn = myXcorr(xn);  

CXk = fft(cxn,nfft);    myCXk = myfft(mycxn,nfft);
Pxx = abs(CXk);         myPxx = abs(myCXk);
index = 0:round(nfft/2-1);
k = index*Fs/nfft;
plot_Pxx = 10*log10(Pxx(index+1));     myplot_Pxx = 10*log10(myPxx(index+1));
subplot(211);    plot(k,plot_Pxx);          title('间接法---内置MATLAB函数');
subplot(212);    plot(k,myplot_Pxx);        title('间接法---自编MATLAB函数');

代码中myfft函数见https://blog.csdn.net/xq_520/article/details/102501908

验证结果（部分截图）：
傅里叶变换：

功率谱密度计算结果：