本文档提供了使用MATLAB实现的经典功率谱估计方法的代码，包括周期图法、Welch平均周期图法和MUSIC算法，适用于信号处理中的功率谱分析。
本段落档提供了三种经典的功率谱估计方法的MATLAB代码：直接法、改进后的直接法（包括Bartlett法）以及Welch法。
**1. 直接法**
也称为周期图法，该方法通过将随机序列x的N个观测数据视作能量有限序列，并计算其离散傅立叶变换X。之后取幅值平方并除以N作为真实功率谱估计。
```matlab
clear;
Fs = 1000; % 设置采样频率为1000Hz
n = 0:1/Fs:1; % 创建时间向量，用于生成含噪声的序列xn
xn = cos(3 * pi * n) + randn(size(n)); % 添加高斯白噪声到信号中
window = boxcar(length(xn)); % 使用矩形窗函数
nfft = 1024;
[Pxx, f] = periodogram(xn, window, nfft, Fs); % 计算功率谱密度估计值Pxx和频率向量f
plot(f,Pxx);
```
**2. 改进的直接法**
对于原始周期图方法，当数据长度N过大时会导致频谱曲线波动增加；而过小则会降低分辨率。改进的方法包括Bartlett平均周期图以及Welch法。
- **Bartlett 法**
Bartlett 法通过将 N 点序列分为若干段计算各自的周期图，并求这些结果的均值，以减少方差。
```matlab
clear;
Fs = 1000; % 设置采样频率为1000Hz
n = 0:1/Fs:1; % 创建时间向量，用于生成含噪声的序列xn
xn = cos(3 * pi * n) + randn(size(n)); % 添加高斯白噪声到信号中
window = boxcar(length(xn));
nfft=1024;
[Pxx, Pxxc] = psd(xn, window, Fs, Fs,NFFT, nfft);
index = 0:round((length(Pxx)-1)/3);
k=index*Fs/nfft; % 计算频率索引
plot_Pxx=10*log10(abs(Pxx)); % 转换为dB值
plot_Pxxc=10*log10(abs(Pxxc));
figure;
plot(f, plot_Pxx);
pause;
figure;
plot(k, plot_Pxxc(index+1));
```
- **Welch 法**
Welch 方法在 Bartlett 方法的基础上进行了两方面的改进：选择适当的窗函数，并允许各段间有重叠，以降低方差。
```matlab
clear;
Fs = 1000;
n=0:1/Fs:1;
xn=cos(3 * pi * n)+randn(size(n)); % 添加高斯白噪声到信号中
window=boxcar(length(xn));
window1=hamming(length(xn)); % 使用汉明窗函数
window2=blackman(length(xn)); % 使用Blackman窗函数
nooverlap = 20;
range=half;
[Pxx,f] = pwelch(xn, window, nooverlap, [], Fs); % 计算功率谱估计值Pxx和频率向量f，使用矩形窗
[Pxx1,f]=pwelch(xn,window1,nooverlap,[],Fs);
[Pxx2,f]=pwelch(xn,window2,nooverlap,[],Fs);
plot_Pxx=10*log10(abs(Pxx)); % 转换为dB值
plot_Pxx1=10*log10(abs(Pxx1));
plot_Pxx2=10*log10(abs(Pxx2));
figure;
plot(f, plot_Pxx);
pause;
figure;
plot(f, plot_Pxx1);
pause;
figure;
plot(f, plot_Pxx2);
```
5星
