MATLAB Answers

Different output from pwelch and my own implementation of welch

16 views (last 30 days)
Jan Slavotinek
Jan Slavotinek on 3 Feb 2017
Answered: Jeremy on 10 Feb 2017
Hi!
I am currently trying to understand and implement how to estimate power spectrum as i need it for further analysis. So far, for analysing of my data i have used pwelch function, which make a nice output. Now i need to implement somehow my own estimate using only fft. So i have studied how welch method works (divide data into more segments (possibly overlaping, using hanning window, and do fft for each and then average the result).
So i have written my own function where i try to do estimate. Problem is, thht output is different from what produce pwelch function.
My welch:
function [freq, ampl] = func_FFT2_welch(s, Fs, parts, OverlapFactor)
FrameLength = 2^(nextpow2(length(s)/parts) -1);
Overlap = int32(OverlapFactor * FrameLength);
Shift = FrameLength - Overlap;
part = s(1 : FrameLength);
s(1: Shift) = [];
tmpFFT = abs(fft(part .* hanning(FrameLength)))/FrameLength;
ampl = 2*tmpFFT(1 : int32(FrameLength/2+1));
cnt = 1;
while(length(s) >= FrameLength)
part = s(1 : FrameLength);
s(1: Shift) = [];
tmpFFT = abs(fft(part .* hanning(FrameLength)))/FrameLength;
ampl = ampl + 2*tmpFFT(1 : int32(FrameLength/2+1));
cnt = cnt +1;
end
ampl = ampl / cnt;
ampl = ampl';
freq = Fs/2*linspace(0,1, int32(FrameLength/2+1));
end
See plot i have enclosed. Blue plot is pwelch output (desired). Red is irrelevant to my question. Yellow is simple fft over whole data array (it has a lot of noise, shape is fuzzy) Purple is output of my implementation of welch (whole data array divided into n frames = 2^(nextpow2(length(s)/parts) -1) long, each windowed and overlap 30%). Noise cancelation is good, but problem is, that the overall shape of spectrum is changed.
for the future use i will need to find a linear approximation of some parts of the spectrum. I dont like the way how my function changed the slope of the spectrum.
Note that result is in dB ... Also note that in this case i have called my function with parts = 50 and OverlapFactor = 0.3;
i call pwelch with same overlap and with various "frame length". This only changes variance, but a slope of a the spectrum is unchanged by this.

  0 Comments

Sign in to comment.

Answers (1)

Jeremy
Jeremy on 10 Feb 2017
At first glance the blue line looks like the square of the purple line. Looking at your code confirms this, pwelch is calculating a PSD which is based on each FFT result x the conjugate of the FFT result so it is a quantity squared result.

  0 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!