function [B,f,t] = dsp_specgram(x,nfft,Fs,win,nover)
% [B,f,t] = dsp_specgram(x,nfft,Fs,win,nover)
% Calculates spectrogram of signal x
% Inputs:
% x     -> x(n)
% nfft  -> number of points for fft calculation
%          nfft > length(win)
% Fs    -> sampling frequency (Hz)
% win   -> scalar for Hanning window length
%          or vector of desired window input
% nover -> number of samples to overlap
%          nover < length(win)
% Outputs:
% B     -> fft matrix [m,n]=size(B)
%          each column j represents
%          the corresponding nfft/2 + 1 point fft for
%          samples at x(a:b)
%          a = (j-1)*(nwin-nover) + 1
%          b = j*nwin - (j-1)*nover
% f     -> coresponding frequency vector for fft rows in col.
%          f = [0:nfft/2]*Fs/nfft;
% t     -> coresponding time vector for each fft column
%          t = [0:n-1]*(length(win)-nover)/Fs
% dsp_specgram.m
% author: Nicholas Kottenstette
% date:   11/9/2004
  
  [mx,mn]=size(x);
  if length(win) > 1
    nwin = length(win);
  else
    nwin = win;
    if nwin == 1
      nwin = nwin+1;
    end
    %Hanning window
    win = (1/2)*(1-cos(2*pi*[0:nwin-1]/(nwin-1)));
  end
  %place x in single column format
  if mx < mn
    x = x';
  end
  if mn > 1
    x = x(:,1);
  end
  if nwin > nfft
    nwin = nfft;
  end
  if nwin <= nover
    nover = nwin/2;
  end
  nfft = 2^ceil(log2(nfft));
  m = (nfft/2)+1;
  n = fix((length(x)-nover)/(nwin-nover));
  B = zeros(m,n);
  f = [0:nfft/2]*Fs/nfft;
  t = [0:n-1]*(length(win)-nover)/Fs;
  xw = zeros(nfft,1);
  %Make sure win is in column format
  [wm,wn]=size(win);
  if wn > wm
    win = win';
  end
  temp = zeros(nfft,1);
  try,
    fp = load('fft_coef/coef.mat');
    c = fp.c;
    fp = load('fft_coef/rind.mat');
    rind = fp.rit{floor(log2(nfft))};
  catch,
    c = [];
    rind = [];
    display('dtmf_specgram.m Warning unable to load coeficients and reverse index');
  end
  for j = 1:n
    a = (j-1)*(nwin-nover) + 1;
    b = a + nwin - 1;
    xw(1:nwin) = win.*x(a:b);
    temp = dsp_fft(xw,nfft,c,rind);
    B(:,j) = temp(1:(nfft/2)+1);
  end