function [aad_spec, aad_P, ad_Q, spec_w_phase] = extract_spec(x,T_wind,T_shift,fs,iFFTSize)
%EXTRACT_SPEC Extract spectrogram of input.
%   [aad_spec, aad_P, ad_Q, spec_w_phase] = extract_spec(x,T_wind,T_shift,fs,iFFTSize)
%
%   Really, a less functional version of the built-in SPECTROGRAM
%   function... but at least I know how it works.
%
%   aad_spec:     magnitude spectrogram
%   aad_P:        time-frequency power
%   ad_Q:         time-frequency power warped with (.)^(1/15) power-law
%   spec_w_phase: complex spectrogram
%
%   X can be a string containing the filename of .wav or .sph file or can
%   be the data itself. T_wind is the window length in seconds, T_shift is
%   the hop size in seconds, fs is the sampling rate in Hz, and iFFTSize is
%   the DFT length in points. The default values are:
%
%   T_wind   = 0.0256; 
%   T_shift  = 0.010;  
%   fs       = 16000;  
%   iFFTSize = 2048;   
%   
%   Written by Mark Harvilla, Spring 2011
%   mharvill@cs.cmu.edu

if ischar(x) == 1
    fid     = fopen(x,'rb');
    if strcmp(x(end-2:end),'wav') == 1 %decide between .wav and .sph
        fseek(fid, 44, 'bof');
    elseif strcmp(x(end-2:end),'sph') == 1
        fseek(fid, 1024, 'bof');
    end
    ad_x = fread(fid,'int16'); %read data, assumes 16-bit uncompressed
    fclose(fid);
else
    ad_x = x;
end

if nargin == 1
    T_wind   = 0.0256; %window width
    T_shift  = 0.010;  %hop size
    fs       = 16000;  %sampling frequency
    iFFTSize = 2048;   %DFT length
end
    
N1 = floor(T_wind*fs);
N2 = floor(T_shift*fs);

w = hamming(N1);                                 %use a hamming window
iNumWindows = 1 + floor((length(ad_x) - N1)/N2); %determine number of windows
aad_spec = zeros((iFFTSize/2) + 1, iNumWindows);
spec_w_phase = aad_spec;

for i_FI = 1:iNumWindows
    x = ad_x( 1 + (i_FI-1)*N2 : N1 + (i_FI-1)*N2 ).* w;
    X = fft(x,iFFTSize); 
    X = X(1:(iFFTSize/2) + 1); %chop redundant part
    
    spec_w_phase(:,i_FI) = X;
    aad_spec(:,i_FI) = abs(X);
end

aad_P = aad_spec .^ 2;  %power
ad_Q = aad_P .^ (1/15); %apply power law nonlinearity to the power