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:
However, when I tried overlapping window like this
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.