Manual implementation of spectogram function

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?

Answer 1

[Pawel Krzyzak]

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.

Answer 2

[Mathieu NOE]

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

[Pawel Krzyzak]

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.