mathieu@14: function [y,e] = powspec(x, sr, wintime, steptime, dither)
mathieu@14: %[y,e] = powspec(x, sr, wintime, steptime, sumlin, dither)
mathieu@14: %
mathieu@14: % compute the powerspectrum and frame energy of the input signal.
mathieu@14: % basically outputs a power spectrogram
mathieu@14: %
mathieu@14: % each column represents a power spectrum for a given frame
mathieu@14: % each row represents a frequency
mathieu@14: %
mathieu@14: % default values:
mathieu@14: % sr = 8000Hz
mathieu@14: % wintime = 25ms (200 samps)
mathieu@14: % steptime = 10ms (80 samps)
mathieu@14: % which means use 256 point fft
mathieu@14: % hamming window
mathieu@14: %
mathieu@14: % $Header: /Users/dpwe/matlab/rastamat/RCS/powspec.m,v 1.3 2012/09/03 14:02:01 dpwe Exp dpwe $
mathieu@14: 
mathieu@14: % for sr = 8000
mathieu@14: %NFFT = 256;
mathieu@14: %NOVERLAP = 120;
mathieu@14: %SAMPRATE = 8000;
mathieu@14: %WINDOW = hamming(200);
mathieu@14: 
mathieu@14: if nargin < 2
mathieu@14:   sr = 8000;
mathieu@14: end
mathieu@14: if nargin < 3
mathieu@14:   wintime = 0.025;
mathieu@14: end
mathieu@14: if nargin < 4
mathieu@14:   steptime = 0.010;
mathieu@14: end
mathieu@14: if nargin < 5
mathieu@14:   dither = 1;
mathieu@14: end
mathieu@14: 
mathieu@14: winpts = round(wintime*sr);
mathieu@14: steppts = round(steptime*sr);
mathieu@14: 
mathieu@14: NFFT = 2^(ceil(log(winpts)/log(2)));
mathieu@14: %WINDOW = hamming(winpts);
mathieu@14: %WINDOW = [0,hanning(winpts)'];
mathieu@14: WINDOW = [hanning(winpts)'];
mathieu@14: % hanning gives much less noisy sidelobes
mathieu@14: NOVERLAP = winpts - steppts;
mathieu@14: SAMPRATE = sr;
mathieu@14: 
mathieu@14: % Values coming out of rasta treat samples as integers, 
mathieu@14: % not range -1..1, hence scale up here to match (approx)
mathieu@14: y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2;
mathieu@14: 
mathieu@14: % imagine we had random dither that had a variance of 1 sample 
mathieu@14: % step and a white spectrum.  That's like (in expectation, anyway)
mathieu@14: % adding a constant value to every bin (to avoid digital zero)
mathieu@14: if (dither)
mathieu@14:   y = y + winpts;
mathieu@14: end
mathieu@14: % ignoring the hamming window, total power would be = #pts
mathieu@14: % I think this doesn't quite make sense, but it's what rasta/powspec.c does
mathieu@14: 
mathieu@14: % that's all she wrote
mathieu@14: 
mathieu@14: % 2012-09-03 Calculate log energy - after windowing, by parseval
mathieu@14: e = log(sum(y));