# 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.

• Welcome to DSP.SE! Can you please include some pictures of why you think it's wrong? e.g. one plot showing the difference between matlab's implementation and yours on an example signal? – Peter K. May 13 '16 at 12:13

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