1function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast)
2% function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast)
3% compute binary spectrum using threshold spectrum as pivot
4% bspectrum = binary spectrum (binary)
5% ps=current power spectrum (float)
6% threshold=threshold spectrum (float)
7% bandfirst = first band considered
8% bandlast = last band considered
9  
10% initialization stuff
11  if( length(ps)<bandlast | bandlast>32 | length(ps)~=length(threshold)) 
12  error('BinDelayEst:spectrum:invalid','Dimensionality error');
13end
14
15% get current binary spectrum
16diff = ps - threshold;
17bspectrum=uint32(0);
18for(i=bandfirst:bandlast)
19  if( diff(i)>0 ) 
20    bspectrum = bitset(bspectrum,i);
21  end
22end
23