Performing PSD from FFT. How to unpack FFT?

Hi everybody,

I would like to ask for help in Signal Processing.

I am trying to perform PSD from FFT, without pwelch, periodogram, etc. because I use my own window function. Window function comprises of 1024 weights, and overlapping that I need to use is 61.9%. My signal has fs*time = 48000*9.23=443040 samples.

I am willing to split my signal into windows of N=1024 samples and weight them with the window function.

    winLength = 2^10;                 %1024
    winTime = winLength / Fs;
    nofWins = floor(noSamples/winLength);
    noWins = (100/38.1)*nofWins - 1;  %this is a gap performing the overlap of 61.9%

    odd = true;

    data = zeros(noWins, 1);
    windowed_data = zeros(noWins, 1);

        for k = 1:(noWins-2)
            j = floor(0.381*k);
            at = j*winLength + 1;
            range = at:(at + winLength - 1);

            if   odd
                 data = y(range, 1);
                 windowed_data = data*mywindow;
                 WD = diag(windowed_data); 
            else
                 data = y(range+winLength, 1);
                 windowed_data = data*mywindow;
                 WD = diag(windowed_data); 
            end;

       odd=~odd;

        %
        % Compute FFT
        %

        spectrum = fft(WD);
        PSD(:, k) = spectrum.^2/N;

end

I am following a book to compute PSD, and the next step after taking FFT there is : “Unpack the result, compute the square magnitude, and average it (separately for each frequency bin).”

I do not understand how can I complete the thing they say. Unpack the result apparently means get amplitudes and frequencies. How?

Compute the square magnitude – of what?

And how can I pack it up together afterwards?

It drives me crazy… Thank you so much for any kind of help, guys!

Good luck to you all,

Elena Che

Answer 1

[Wayne King]

Hi Elena, You can use your own window function with pwelch.m.

You just supply a vector, then you can just overlap by 634 samples for the NOVERLAP argument that will give you 61.91% overlap

   x = randn(443040,1);
   win = hamming(1024);
   [Pxx,Fxx] = pwelch(x,win,634);

You can of course specify additional input arguments like NFFT and the sampling rate. For example, assuming a sampling rate of 1 kHz.

   t = 0:0.001:(443040*0.001)-0.001;
   x = cos(2*pi*100*t)+randn(size(t));
   [Pxx,Fxx] = pwelch(x,window,634,length(window),1000);
   plot(Fxx,10*log10(Pxx)); xlabel('Hz'); ylabel('dB');

[Elena]

Hi Wayne,

thanks for your answer. Actually, and unfortunately I cannot use my window function in pwelch.m. At first, because as I remember, pwelch has a Hamming window as a basis, I, in turn, have a window function based on Kaiser, adjusted in the way to increase the peak-sidelobe level (PSLL); at second, my window function is perfectly adjusted to my noise I am measuring (in dB); at third, because I am also going to delete DC average in the noise; and at last... because my professor WANTS me to fo so. AH! :)

[Wayne King]

Hi Elena, pwelch() does not use a Hamming window if you supply your own window. It uses a Hamming window by default, but if you specify a window that overrides that and uses your window. You can use buffer() to divide your signal up into overlapping vectors (as columns of a matrix) and do it that way.