%%% BLIND INVERSE GAMMA CORRECTION%%%    1. generate random signal from N(0,1) of length 'dim'%%%    2. rescale signal into [0,1]%%%    3. apply gamma correction (g(u) = u^gamma)%%%    4. blindly estimate gamma by minimizing bicoherence%%%       for each (possibly overlapping (wininc)) windowed signal (windim)%%%          for each inverse gamma value in 'range'%%%             apply inverse gamma %%%             compute bicoherence of inverse gamma signal%%%          end%%%          estimated gamma = inverse gamma with minimum bicoherence%%%       end%%%    5. final estimate = average of individual estimates%%%%%% RUN:%%%     gammaest = blindgamma_orginal( 1.2, [0.5:0.1:2.0] );%%%%%% REFERENCES:%%%    www.cs.dartmouth.edu/farid/publications/ip01.html%%%    www.cs.dartmouth.edu/farid/publications/josa01.html%%%function[gammaest] = blindgamma( gamma, range )dim    = 2^13;               % signal lengthwindim = 2^10;               % window lengthwininc = 2^9;                % window increment%randn('state',1);f      = randn( 1, dim );    % original signal%f = rand(1, dim);%f = randperm(dim);f      = f - min(f);         % normalize into [0,1]f      = f / max(f);         % fg     = f.^gamma;           % apply gamma %fg = f;ind = [1:windim];k = 1;while( max(ind)<dim )        % for each window   fw = fg(ind);                 for g = 1 : length(range) % for each inverse gamma value      bic = blindgamma( fw.^(1/range(g)) );      B(k,g) = mean( abs(bic(:)) );   end   k = k + 1;   ind = ind + wininc;endfor k = 1 : size(B,1)   [val,ind] = min(B(k,:));   gammaest(k)  = range(ind);endgammaest = mean(gammaest);return;%%%%%% COMPUTE BICOHERENCE%%%function [bic] = bispec( y )[ly, nrecs] = size(y);if( ly == 1 )   y = y(:);end% ---------------- bispectrum parameters ----------------------nsamp   = 64;overlap = 32;nadvance = nsamp - overlap;nrecs   = fix ( (ly*nrecs - overlap) / nadvance);nfft    = 128;wind    = hanning(nsamp); wind    = wind(:);% ---------------- bispectrum calculation ----------------------bic	= zeros( nfft, nfft );YY 	= zeros( nfft, 1 );Yc12	= zeros( nfft, nfft );YY12	= zeros( nfft, nfft );mask	= hankel([1:nfft],[nfft,1:nfft-1] );ind	= [1:nsamp];for k = 1 : nrecs        ys	= y(ind);	ys	= (ys(:) - mean(ys)) .* wind;	ys 	= ys(:) - mean(ys);	Y	= fft(ys,nfft) / nsamp;        Yc	= conj( Y );	YY	= YY + ( Y .* Yc );	Yc12(:) = Yc( mask );	bic	= bic + (Y * Y.') .* Yc12;        ind	= ind + nadvance;endYY      = sqrt( YY / nrecs );YY12(:) = YY(mask);bic     = bic / nrecs;bic     = bic ./ sqrt( (abs(YY * YY.').^2) .* (abs(YY12).^2) );bic     = fftshift( bic );return;