why does spectrogram of stft shows different magnitude while changing the window size?
21 views (last 30 days)
Show older comments
fs=5000; % Sampling rate
N=1000; % Number of data points
y=sin(500*pi*[0:1:N-1]/fs)...
+2*sin(500*1.5*pi*[0:1:N-1]/fs)...
+3*sin(500*3*pi*[0:1:N-1]/fs)...
+4*sin(500*4*pi*[0:1:N-1]/fs);t=[0:1:N-1]/fs;
figure(1);
windowsize = 200;
window = boxcar(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')
For the above signal, the expected amplitudes of frequencies are in the range 1-4. once we i do the stft with window size 200, i am getting amplitudes of frequency in the range of 0-400 Hz. if i reduce the window size to 16, the amplitudes moving to the range of 0-45 Hz. Did i do anything wrong in my code or am i conceptually wrong somewhere?
0 Comments
Accepted Answer
Mathieu NOE
on 17 Mar 2021
hello
I was tired of not being able to understand why some matlab functions (mostly based on fft) would not scale the output as I was expecting.
Finally decided (long time ago) to rewritte the spectrogram function myself - now the amplitude of the spectrogram are ok with the input data. I only decided to simplify the function, so it's only hanning window for the time being. If you change the type of window , you have to correct for the new correction coefficient.
also the overlap is a percentage of nfft (most of the time I leave it at 75%)
hope it helps
fs=5000; % Sampling rate
N=1000; % Number of data points
t=[0:1:N-1]/fs;
y=sin(500*pi*t)...
+2*sin(500*1.5*pi*t)...
+3*sin(500*3*pi*t)...
+4*sin(500*4*pi*t);
figure(1);
windowsize = 200;
window = boxcar(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
% [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); % S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs)
[S,F,T] = myspecgram(y, fs, nfft, 0.75); % overlap is 75% of nfft here
imagesc(T,F,(abs(S)))
colorbar('vert');
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')
colormap('jet');
function [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal (example sinus amplitude 1 = 0 dB after fft).
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,1);
s_tmp((1:samples)) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_specgram = [];
for ci=1:spectnum
start = (ci-1)*offset;
sw = signal((1+start):(start+nfft)).*window;
fft_specgram = [fft_specgram abs(fft(sw))*4/nfft]; % X=fft(x.*hanning(N))*4/N; % hanning only
end
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;
% time vector
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;
end
0 Comments
More Answers (0)
See Also
Categories
Find more on Time-Frequency Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!