信息处理 - 实现 Welch 的功率谱密度方法 - 吾爱随笔录

我想在 MATLAB 中实现韦尔奇的 PSD 计算方法。我不想使用内置的 MATLABcpsd函数，这样我就可以更改 FFT 实现cpsd（出于某种目的）。

我已经尝试根据 Welch 论文和页面上的解释复制该方法：非参数方法。下面是我已经实现的代码，其中包含关于实现结果的步骤的注释。（它用于加速度计的真实数据，所以numaxis代表轴数， $xyz$ . 每列代表各个轴，时间序列数据按行排列）：

function psd = pspectra(data1,data2, window, overlap, nfft, fs)

%checking the size
if(size(data1,2) ~= size(data2,2))
    error('size of column is not equal');
elseif(size(data1,1) ~= size(data2,1))
    error('size of row is not equal');
else
    numaxis = size(data1,2);
end

%calculating window sliding step for iteration
winsize = length(window);
step = winsize - overlap;
iter = 1 + (size(data1,1) - winsize)/step;

%start and end index of first window/segment
istart = 1;
iend = istart + winsize - 1;

%start calculating fft for each window
for i=1:iter
    for j=1:numaxis
        %apply window
        data1(istart:iend,j) = data1(istart:iend,j).*window;
        data2(istart:iend,j) = data2(istart:iend,j).*window;

        %calculate fft
        fft1(:,j,i) = fft(data1(istart:iend,j),nfft);
        fft2(:,j,i) = fft(data2(istart:iend,j),nfft);     
    end

    %move to next window segment
    istart = istart + step;
    iend = iend + step;
end

%obtain scale to create modified periodogram
scale = 1/(fs.*sum(window.*window));

%averaging window result and apply the scaling
psd = zeros((nfft/2)+1,size(fft1,2));

for i=1:iter
    psd = psd + fft1(1:(nfft/2)+1,:,i).*conj(fft2(1:(nfft/2)+1,:,i));
end

psd = psd.*scale./iter;

%multiply by 2 except dc and nyquist component (1 and 51)
psd(2:50,:) = psd(2:50,:).*2;

end

cpsd()我尝试使用具有非重叠窗口/段的多个段，并得到与 MATLAB 函数完全相同的结果。例如，使用此调用：

pspectra(data1,data1,hamming(50),0,100,100);
pspectra(data1,data1,hamming(25),0,100,100);

但是，当我尝试像这样重叠窗口时

pspectra(data1,data1,hamming(50),25,100,100);
pspectra(data1,data1,hamming(60),20,100,100);

得到的功率谱完全不同。我试图找出我的实现是否错误（理论上或代码），不幸的是还没有运气。我想念这里的东西吗？

我真的很感激任何帮助。