MATLAB Answers

Manual implementation of spectogram function

14 views (last 30 days)
Hello,
I am trying to implement my own function that gives the same results as Matlab spectogram function.
So far I have accomplished function like this:
function out = manulaSpectogram(x, win, noverlap, nfft)
x = x(:);
n = length(x);
wlen = length(win);
nUnique = ceil((1+nfft)/2); % number of uniqure points
L = fix((n-noverlap)/(wlen-noverlap)); % number of signal frames
out = zeros(L, nUnique);
index = 1:wlen;
for i = 0:L-1
xw = win.*x(index);
X = fft(xw, nfft);
out(i+1, :) = X(1:nUnique);
index = index + (wlen - noverlap);
end
end
In my tests it works perfectly and gives the same results like spectogram function when parameter nfft is greater or equal to length of window.
% first test (nnft = window length):
A = [1,2,3,4,5,6];
window = 6;
overlap = 2;
nfft = 6;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
9.7300 + 0.0000i -5.2936 + 0.9205i 0.7279 - 0.3737i -0.1186 + 0.0000i
s2 =
9.7300 + 0.0000i -5.2936 - 0.9205i 0.7279 + 0.3737i -0.1186 + 0.0000i
% second test (nfft > window length):
A = [1,2,3,4,5,6];
window = 3;
overlap = 2;
nfft = 6;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
2.3200 + 0.0000i 0.9600 + 1.9399i -1.0400 + 1.5242i -1.6800 + 0.0000i
3.4800 + 0.0000i 1.5000 + 2.8752i -1.5000 + 2.3209i -2.5200 + 0.0000i
4.6400 + 0.0000i 2.0400 + 3.8105i -1.9600 + 3.1177i -3.3600 + 0.0000i
5.8000 + 0.0000i 2.5800 + 4.7458i -2.4200 + 3.9144i -4.2000 + 0.0000i
s2 =
2.3200 + 0.0000i 0.9600 - 1.9399i -1.0400 - 1.5242i -1.6800 + 0.0000i
3.4800 + 0.0000i 1.5000 - 2.8752i -1.5000 - 2.3209i -2.5200 + 0.0000i
4.6400 + 0.0000i 2.0400 - 3.8105i -1.9600 - 3.1177i -3.3600 + 0.0000i
5.8000 + 0.0000i 2.5800 - 4.7458i -2.4200 - 3.9144i -4.2000 + 0.0000i
In the case when length of window is less than nfft than the results are different.
% third test (nfft < window length):
A = [1,2,3,4,5,6];
window = 6;
overlap = 2;
nfft = 3;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
9.7300 + 0.0000i 0.7279 - 0.3737i
s2 =
3.6121 + 0.0000i -1.6861 + 1.6807i
So how can I improve my function to recieve the same results even in the case when nnft is less than window length? How Matlab's spectogram calculates this case?

Accepted Answer

Pawel Krzyzak
Pawel Krzyzak on 29 Apr 2021
Edited: Pawel Krzyzak on 29 Apr 2021
I noticed that when window size is greater than nfft scalar number, the data has to be transformed somehow. Finally I found an inner Matlab function that probably is called in the original spectogram Matlab function. It is named datawrap and wraps input data modulo nfft.
So in my function I had to transform data segment (in the same way how datawrap function does it) before calling fft. Improved function:
function out = manulaSpectogram(x, win, noverlap, nfft)
x = x(:);
n = length(x);
wlen = length(win);
nUnique = ceil((1+nfft)/2); % number of uniqure points
L = fix((n-noverlap)/(wlen-noverlap)); % number of signal frames
out = zeros(L, nUnique);
index = 1:wlen;
for i = 0:L-1
xw = win.*x(index);
% added transformation
if length(xw) > nfft
xw = sum(buffer(xw, nfft), 2);
end
% end of added transformation
X = fft(xw, nfft);
out(i+1, :) = X(1:nUnique);
index = index + (wlen - noverlap);
end
end
I believe it works properly because it gives the same results as Matlab spectogram function.

More Answers (1)

Mathieu NOE
Mathieu NOE on 29 Apr 2021
hello
I share with you my own version , hope it helps !
notice i always use hanning window and the correction factor for amplitude is related to this choice
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
  1 Comment
Pawel Krzyzak
Pawel Krzyzak on 29 Apr 2021
Thank your very much for your implementation but I still don't know how can I modify my function to recieve the same results as results given by matlab spectogram funcion.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!