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.