Re-implementing the “spectrogram” function from matlab

I am trying to make a spectrogram viewer without using the spectrogram command. For this purpose, I used a sin function. I broke up the input signal into 256 segments, multiplied each segment with the Hamming window, took the FFT of each slice and then stored them to obtain spectra of the signal.

My problem is that I couldn't get the exact spectrum. It is probably due to the fact that I could not divide input signal into overlapping segments?

This is my code: Can someone help me where I am doing it wrong?

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 10240;                     % Length of signal
NFFT=10240;                   % number of fft points
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = hann(L).*x';
figure
plot(X)
k=256;

for idx = 1:length(x)/k
%take a slice from the input vector
slice =  x(1+(idx-1)*k:idx*k);
%multiply it with the window and transform it into frequency domain
spectrum =  fft(slice'.*hamming(k),k);
%get the spectrum magnitude at each of the 256 frequency points and store it
mag_spectrum(:,idx*1024) = abs(spectrum).^2;

end

figure
w=hamming(k);
plot(w)
figure
plot(t,X,'k');
xlabel('Time (sec)','FontSize',14);
ylabel('\it{x(t)}','FontSize',14);

% spectrogram(x,128);
% SET PARAMETERS

figure
plot(slice); % To look for it is windowed or not which is not in my case?

%figure
%plot(mag_spectrum);
figure
spectrogram(X,256); % Real spectogram


Differences between your implementation and what the spectogram function is displaying.

1. 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'

2. You don't have the proper scaling (see reference paper and code below). Also, I believe the spectogram is displayed in dB.

3. 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 < 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) > 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

• I suppose you did not understand my problem excatly I am not looking for a new code. I just want to know where I am doing it wrong? Thank you anyway – Serkan Dec 8 '13 at 23:16
• Can you elaborate on what you mean by "My problem is that I couldn't get the exact spectrum." What is the 'exact spectrum' supposed to look like? – porten Dec 9 '13 at 2:35
• Spectogtam of the signal. I compare my result with spectogram but they are different. – Serkan Dec 9 '13 at 9:38
• Thank you I make some modificaitons on my code with your opinions and it is better now ;) – Serkan Dec 9 '13 at 19:16
• The link that you provide did not work paper for scaling. I plot spectogram but my time and frequencies are normalized I cant see the where real frequency value.I try to multiply mag_spectrum values but It did not work – Serkan Dec 19 '13 at 20:10