comparison nonExposed/utils/powspec.m @ 14:b1901e8d8f5f

initial commit
author Mathieu Lagrange <mathieu.lagrange@cnrs.fr>
date Tue, 17 Mar 2015 09:34:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 14:b1901e8d8f5f
1 function [y,e] = powspec(x, sr, wintime, steptime, dither)
2 %[y,e] = powspec(x, sr, wintime, steptime, sumlin, dither)
3 %
4 % compute the powerspectrum and frame energy of the input signal.
5 % basically outputs a power spectrogram
6 %
7 % each column represents a power spectrum for a given frame
8 % each row represents a frequency
9 %
10 % default values:
11 % sr = 8000Hz
12 % wintime = 25ms (200 samps)
13 % steptime = 10ms (80 samps)
14 % which means use 256 point fft
15 % hamming window
16 %
17 % $Header: /Users/dpwe/matlab/rastamat/RCS/powspec.m,v 1.3 2012/09/03 14:02:01 dpwe Exp dpwe $
18
19 % for sr = 8000
20 %NFFT = 256;
21 %NOVERLAP = 120;
22 %SAMPRATE = 8000;
23 %WINDOW = hamming(200);
24
25 if nargin < 2
26 sr = 8000;
27 end
28 if nargin < 3
29 wintime = 0.025;
30 end
31 if nargin < 4
32 steptime = 0.010;
33 end
34 if nargin < 5
35 dither = 1;
36 end
37
38 winpts = round(wintime*sr);
39 steppts = round(steptime*sr);
40
41 NFFT = 2^(ceil(log(winpts)/log(2)));
42 %WINDOW = hamming(winpts);
43 %WINDOW = [0,hanning(winpts)'];
44 WINDOW = [hanning(winpts)'];
45 % hanning gives much less noisy sidelobes
46 NOVERLAP = winpts - steppts;
47 SAMPRATE = sr;
48
49 % Values coming out of rasta treat samples as integers,
50 % not range -1..1, hence scale up here to match (approx)
51 y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2;
52
53 % imagine we had random dither that had a variance of 1 sample
54 % step and a white spectrum. That's like (in expectation, anyway)
55 % adding a constant value to every bin (to avoid digital zero)
56 if (dither)
57 y = y + winpts;
58 end
59 % ignoring the hamming window, total power would be = #pts
60 % I think this doesn't quite make sense, but it's what rasta/powspec.c does
61
62 % that's all she wrote
63
64 % 2012-09-03 Calculate log energy - after windowing, by parseval
65 e = log(sum(y));