Different output from pwelch and my own implementation of welch

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.

Answer 1

[Jeremy]

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.