matlab代码，修正协方差功率谱密度估计算法

MATLAB是一种强大的数值计算环境，常用于信号处理和数据分析。在估计协方差功率谱密度（Covariance Spectral Density，CSD）时，尤其在时间序列分析或频域分析中，可能需要使用傅立叶变换或其他统计方法。一种常见的改进方法是对传统的快速傅立叶变换（FFT）结果进行平滑或窗函数处理，以减少噪声影响。

以下是简单的步骤概述：

% 假设你有一个包含时间序列数据的向量x
data = ...; % 你的实信号数据

% 定义窗口函数，如汉明窗（Hamming window）
window = hamming(length(data));

% 应用窗口到数据上
windowed_data = data .* window;

% 计算窗函数加权的数据的快傅立叶变换
csd_est = abs(fft(windowed_data)) ./ length(windowed_data); % 分母除以长度是为了归一化

% 对谱估计进行平滑，例如使用移动平均
window_length_smooth = 50; % 你可以调整这个参数
csd_smoothed = movmean(csd_est, window_length_smooth);

% 将谱估计转换回频率轴，假设采样率fs
[~, freq] = freqz(ones(1, length(window)), [1, -1], length(data)/length(window)); % 使用fftshift不是必需的，取决于你的需求

% 修正低频部分可能出现的泄漏效应
csd_corrected = csd_smoothed ./ (1 + (freq / (fs/2))^2); % fs/2是Nyquist频率

% 保存或显示结果
figure;
plot(freq, csd_corrected);
xlabel('Frequency (Hz)');
ylabel('Covariance Power');
title('Corrected CSD Estimate');