Special case of spectrogram not match FFT

3 Ansichten (letzte 30 Tage)
John
John am 31 Aug. 2018
Bearbeitet: William Rose am 29 Jan. 2024
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

Antworten (1)

William Rose
William Rose am 29 Jan. 2024
Bearbeitet: William Rose am 29 Jan. 2024
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.
  2 Kommentare
William Rose
William Rose am 29 Jan. 2024
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
William Rose am 29 Jan. 2024
Bearbeitet: William Rose am 29 Jan. 2024
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.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Time-Frequency Analysis finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by