Special case of spectrogram not match FFT

Hi,

I am comparing a special case of spectrogram/short-time Fourier transform with FFT. I was expecting that spectrogram of the signal, which has only one window covering the whole signal duration, will match FFT. The FFT is single-sided spectrum. I could match FFT result except at DC if dividing STFT with its size. At DC, it looks like the result is double of FFT result. The code is below. I have two questions with this:

1) Why do I have to divide STFT with its size? Wouldn't it be that the units of both STFT and FFT are same as unit of the signal in time domain?

2) Why after doing division, it doesn't match for DC case?

Thanks,

if true

 clear all;
 clc;
 close all;
 nfft = 1024;
 Fs = 100000;
 Ts=1/Fs;
 time_vect = 0:Ts:(nfft-1)*Ts;
 signal_vect = sin(time_vect/Ts*0.1)+0.5;
 [S_signal,F_signal_vect,T_vect_signal] = spectrogram(signal_vect,hann(nfft),0,nfft,Fs);
 abs_S_signal = abs(S_signal);
 Y = fft((signal_vect(1:nfft)).*hanning(nfft)');
 f_fft = Fs*(0:(nfft/2))/nfft;
 P2 = abs(Y/nfft);
 freq_content_signal_FFT = P2(1:nfft/2+1);
 freq_content_signal_FFT(2:end-1) = 2*freq_content_signal_FFT(2:end-1);
 figure;
 subplot(121);
 plot(time_vect,signal_vect,'b');
 title('Signal vs time');
 xlabel('Time[s]');
 subplot(122);
 hold on;
 plot(f_fft,freq_content_signal_FFT,'b','DisplayName','FFT');
 plot(F_signal_vect,abs_S_signal/numel(abs_S_signal),'r','DisplayName','SFFT/size');
 legend show;
 title('Frequency content vs. freq');

end

Answer 1

[William Rose]

@John,

There are many ways to normalize a spectrum. Mathematicians like to have a factor 1/sqrt(n) in the forward and reverse transforms, because it makes them symmetric. For a power spectrum, it is customary to divide the DFT (which is what the FFT algorithm computes) by signal length N, because if you don't, then a signal that is twice as long will have spectrum values that are twice as large. If you want the power spectral density, you also divide by the frequency spacing, delta f. The discrepancy you observed at DC is because the one sided power spectrum of a real sequence* is related to the 2 sided by

Pxx1(f)=Pxx2(f) when f=0 and when f=fs/2

Pxx1(f)=2*Pxx2(f) when 0<f<fs/2

*For a complex sequence, the spectrum is generally not symmetric about the Nyquist frequency (or, equivalently, about f=0), so it is not appropriate to compute a 1-sided spectrum.

[William Rose]

@John,

Here is a demo of the relationship between th DFT and the 2-sided and 1-sided power spectra computed with periodogram.

rng(1);             % seed random number generator, for reproducibility
x=randn(8,1);       % random sequence
X=fft(x);           % X=discrete Fourier transform
Pxx2=periodogram(complex(x),ones(8,1),8,'power');   % 2-sided spectrum
Pxx1=periodogram(x,ones(8,1),8,'power');            % 1-sided spectrum
% put results in a table to compare them
T=table(X,(abs(X)/8).^2,Pxx2,[Pxx1;[0;0;0]],'VariableNames',["X(f)","(|X|/N)^2","Pxx2","Pxx1"])

T = 8×4 table
          X(f)           (|X|/N)^2      Pxx2        Pxx1  
    _________________    _________    ________    ________

     -3.1344+0i           0.15351      0.15351     0.15351
      2.3474-0.12962i    0.086363     0.086363     0.17273
    -0.17743-1.5397i     0.037535     0.037535    0.075071
     -1.9544-0.52917i    0.064055     0.064055     0.12811
      -2.489+0i          0.096796     0.096796    0.096796
     -1.9544+0.52917i    0.064055     0.064055           0
    -0.17743+1.5397i     0.037535     0.037535           0
      2.3474+0.12962i    0.086363     0.086363           0

The table above shows that the 2-sided power spectrum Pxx2(f) = (|X(f)|/N).^2. It shows that the one-sided spectrum has double the power for 0<f<fs/2.

[William Rose]

@John,

Spectrogram's spectrum, s, is identical to X(f) returned by fft(x), if one calls spectrogram() with the appropriate inputs:

rng(1);             % seed random number generator, for reproducibility
x=randn(8,1);       % random sequence
X=fft(x);           % X=discrete Fourier transform
[s,~,~]=stft(x,Window=rectwin(8),FFTLength=8,FrequencyRange="twosided");
disp([X,s]);  % compare X, s

  -3.1344 + 0.0000i  -3.1344 + 0.0000i
   2.3474 - 0.1296i   2.3474 - 0.1296i
  -0.1774 - 1.5397i  -0.1774 - 1.5397i
  -1.9544 - 0.5292i  -1.9544 - 0.5292i
  -2.4890 + 0.0000i  -2.4890 + 0.0000i
  -1.9544 + 0.5292i  -1.9544 + 0.5292i
  -0.1774 + 1.5397i  -0.1774 + 1.5397i
   2.3474 + 0.1296i   2.3474 + 0.1296i

They match.

If you use a hann window...

X=fft(x.*hann(8,'periodic'));           % X=discrete Fourier transform
[s,~,~]=stft(x,Window=hann(8,'periodic'),FFTLength=8,FrequencyRange="twosided");
disp([X,s]);  % compare X, s

  -2.7409 + 0.0000i  -2.7409 + 0.0000i
   2.0017 + 0.3201i   2.0017 + 0.3201i
  -0.1870 - 0.6052i  -0.1870 - 0.6052i
  -0.3106 + 0.1203i  -0.3106 + 0.1203i
  -0.2673 + 0.0000i  -0.2673 + 0.0000i
  -0.3106 - 0.1203i  -0.3106 - 0.1203i
  -0.1870 + 0.6052i  -0.1870 + 0.6052i
   2.0017 - 0.3201i   2.0017 - 0.3201i

They are identical in this case also. I wasn't sure how that was going to work out. If you estimate the power spectrum with Matlab's pwelch() and a non-rectangular window, then the equality between X=fft(x.*window) and the spectral estimate is not so perfect, becaue pwelch() does an internal adjustment to make up for the power loss due to windowing.