Special case of spectrogram not match FFT

5 次查看（过去 30 天）
John
John 2018-8-31
编辑： William Rose 2024-1-29
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

请先登录，再回答此问题。

回答（1 个）

William Rose
William Rose 2024-1-29
编辑：William Rose 2024-1-29
@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.
  2 个评论
William Rose
William Rose 2024-1-29
Ran in:
@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
William Rose 2024-1-29
编辑：William Rose 2024-1-29
Ran in:
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by