Re-implementing the "spectrogram" function from matlab
Differences between your implementation and what the spectogram function is displaying.


This looks wrong and you don't preallocate the mag_spectrum matrix (which is bad practice).

mag_spectrum(:,idx*1024) = abs(spectrum).^2;

mag_spectrum(idx, :) = abs(spectrum).^2; % Use this instead

mag_spectrum(:, idx) = abs(spectrum).^2; % Or this... depending on how you want the data to 'waterfall'
You don't have the proper scaling (see reference paper and code below). Also, I believe the spectogram is displayed in dB.
You're saving the two-sided FFT. This is unnecessary and the spectogram defaults to the one-sided FFT. Use the proper scaling accordingly.


Here's some code used to compute the average PSD of an input sequence using a provided window and some other parameters. You can modify the last few lines to return the entire spectogram by getting rid of the call to 'mean' and properly handling the last step of reindexing.

Performing the overlap is pretty easy if you convert the data vector into a matrix using reshape, and then proceed to concatenate that matrix with a shifted version of itself.

A good explanation of some of the subtleties of the PSD is presented in this paper, especially regarding scaling.

[px f] = welch(randn(1,32768), 1024, ones(1,1024), 512, 1024, 0);
plot(f, 10*log10(abs(px)))

function [Px, f] = welch(x, fs, win, overlap, Nfft, oneside)
% x=vector of data samples
% fs=sampling frequency (in Hz)
% win=vector containing the window
% overlap=number of points to overlap
% Nfft=length of FFT to use
% oneside=a logical variable indicating whether 
%    scaling is for 1-sided or 2-sided PSD (true if 1-sided)
% Px=estimate of the power spectral density
% f=vector of frequencies (in Hz) for plotting
if nargin &lt; 6
  error('Insufficient input arguments')
end

x = x(:);
win = win(:);

Nw = length(win);
Ns = length(x);

% Compute the frequency vector and scale factor
if oneside == true
  sfactor = 2 ./ (fs * sum(win.^2));
  f = fs * (0:Nfft/2).'./Nfft;
else
  sfactor = 1 ./ (fs * sum(win.^2));
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end

% Clip excess data
x = x(1:floor(Ns/(Nw-overlap)) * (Nw-overlap));

% Perform overlap by appending a shifted version
% of the data matrix
y = reshape(x, Nw-overlap, []);
y = [y; circshift(y(1:overlap, :), [0 -1])];
y = bsxfun(@times, y, win);

% get rid of last column
if size(y, 2) &gt; 1
  y = y(:,1:end-1);
end

% Compute the PSD
Px = sfactor * mean( abs(fft(y, Nfft)).^2, 2);

% Return full or onesided result
if oneside == true
  Px = Px(1:end/2+1);
else
  Px = fftshift(Px);
end