Compute PSD from fft compared with spectrogram

Hi there,

Recently I am working on translating some Matlab code to Java, as I want to do some pattern recognition in mobile phones. In the Matlab code from my teammate, he used this function "spectrogram" to get PSD, which I couldn't find a complete open-source implementation online. Therefore I have to do it by my own.(who doesn't!!)

I got some algorithms from the forum such as:

y = data;
Fs = 90;
L = length(y);
y  = y(:) - mean(y);
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT);
f = Fs/2*linspace(0,1,NFFT/2+1);
psd = abs(Y(1:NFFT/2+1)).^2/L;
psd = psd/Fs;

which has the same result with

   periodogram

But NONE of them has the same result as the spectrogram. Would any one help me out here? Thank you very much.

Answer 1

[Wayne King]

spectrogram optionally returns the short-time periodogram estimates. I'll show you how to establish the correspondence for the default Hamming window used in spectrogram. I'll establish the correspondence only for the 1st segment in the following example, but it is straightforward to extend it over the entire data record with the proper overlap.I'll just use a noise vector and assume a sampling rate of 1 kHz.

   x = randn(1000,1);
   [y,f,t,p] = spectrogram(x,256,250,256,1000);
   % now to create an equivalence between the column vectors
   % in p
   win = hamming(256);
   xdft = fft(x(1:256).*win);
   psdx = abs(xdft(1:length(xdft)/2+1)).^2;
   psdx = 1/(1000*norm(win,2)^2)*psdx;
   psdx(2:end-1) = 2*psdx(2:end-1);
   % now compare
   subplot(211)
   plot(psdx)
   subplot(212)
   plot(p(:,1))

[Yifei]

Hi Wayne,

Sorry that I was on my final exam and didn't have chance to look into your code until now. Thank you so much for your help. It is amazing. The code is so neat and makes perfect sense to me.

However as I noticed, if the input parameter for F is a vector, the "spectrogram" function will compute at the frequency F with Goertzel Algorithm. And it becomes one-sided. What's the difference here in algorithm sense?

Anyway, thanks to your help, I have figures out how to solve my problem. And I also learnt a lot during this process.

[Brandon Madsen]

This is a very helpful explanation, but I am unsure about the purpose of a certain line of code, namely: psdx = 1/(1000*norm(win,2)^2)*psdx; . Could someone explain what this part is doing and why it is needed?

EDIT: Oh, nm, it's for normalization to Hz I assume.

[HYZ]

If I want to plot for different frequencies, let's say 0 to 300Hz in the y-axis [y,f,t,p] = spectrogram(x,256,250,1:300,1000); how can I change the below?

win = hamming(256);

xdft = fft(x(1:256).*win);

psdx = abs(xdft(1:length(xdft)/2+1)).^2;

psdx = 1/(1000*norm(win,2)^2)*psdx;

psdx(2:end-1) = 2*psdx(2:end-1);

Answer 2

[Vieniava]

Spectrogram is time-frequency (3D=time vs freq. vs amplitude) representation of a signal and periodogram/fft is frequency only (2D= freq vs amplitude) representation. Spectrogram shows how the frequency spectrum is changing over the time. Spectrogram is a set of consecutive fft's. Spectrogram is a matrix and fft/peridogram is a vector.

[Yifei]

Hi,

Thanks for the answer first.

Sorry for the vague description about my question. Actually what I am looking for is the mathematical implementation for spectrogram PSD via standard fft.

After looking into documents and Matlab implementations, I found that PSD in spectrogram is x^2 of the STFT, which is based on DFT via Goertzel algorithm over hamming window. Honestly I am still a little confused about the algorithm and will continue working on that later(have one final tmrw, so...). If you have ideas about how and why it is computed this way, please let me know.

Thank you for your time~