matlab - Implementing Welch's method for Power Spectral Density

I want to implement Welch's method for PSD calculation in MATLAB. I do not want to use built in MATLAB cpsd function, such that I can change the FFT implementation in cpsd(for certain purpose).

I already tried to replicate the method based on Welch paper and explanation from the page: Nonparametric Methods. Below is the code that I already implemented with comments about steps to achieve the result. (it is used for accelerometer real data, so numaxis represent number of axis, $xyz$. Each column represent each axis, the time series data is arranged into row):

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

I have tried using multiple segment with non overlapping window/segment, and got exactly the same result as cpsd() function from MATLAB. For example, using this call:

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

However, when I tried overlapping window like this

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

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

The resulted power spectra is totally different. I have tried to find if my implementation is wrong (in theory or code), and unfortunately no luck yet. Do I miss something in here?

I really appreciate for any help.

Answer

I believe these lines are you problem:

    data1(istart:iend,j) = data1(istart:iend,j).*window;
    data2(istart:iend,j) = data2(istart:iend,j).*window;

you are modifying the incoming data. When you don't have overlapping windows, you don't notice it. However, when you have overlapping windows, the next window will include data that has already had a window applied... which will not be what you want.