clear all;
figure(1);clf;
%load handel; fs=8000;
%[y,fs]=wavread('C:\Documents and Settings\Bruce\My Documents\Koe\chickadee.wav');
%[y,fs]=wavread('C:\Documents and Settings\Bruce\My Documents\Koe\voicesample.wav');
[y,fs]=wavread('C:\Documents and Settings\Bruce\My Documents\Koe\woodthrush.wav');
y = y(1:fs*4);
[b,a] = butter(6,.01,'high')
y = filter(b, a, y);

fs = 10000;
 t=0:1/fs:1;                 % 2 secs @ 1kHz sample rate
 y = chirp(t,0,1,4000);
 %y= sin(2*pi*625*t);
y = -1 + 2*(y>0);

figure(1);clf;
set(gcf,'doublebuffer','on')

%minimum p = 4
p = 7 ;
if p<4
    error('requires >= 16 points')
end
n = 2^p ;
%generate sequency order
% Sn(2k)=Sn-1(k) and Sn(2k+1)=(2^n)-1-Sn-1(k)
s{3} = [0 7 3 4 1 6 2 5];
for i=4:p
    s{i} = zeros(1,2*length(s{i-1}));
    s{i}(1:2:2^i) = s{i-1};
    s{i}(2:2:2^i) = 2^i-1 - s{i-1};
end

%loop thru waveform
slice=1;
for tt=1:n:length(y)-n
    x = y(tt:tt+n-1);% .* hanning(n)' ;
    %natural order FWT
    bb = n/2; st = n;
    for i=1:p
        for j=1:st:n
            for k=0:bb-1
                ii=j+k;
                t1 = x(ii); t2 = x(ii+bb);
                x(ii) = t1+t2; x(ii+bb) = t1-t2;
            end
        end
        bb=bb/2; st=st/2;
    end
    %normalize
    x = x/n;
    %and reorder the spectrum
    x(s{p}+1) = x;
    %and combine sal and cal components 
    %first index is DC
    %last is highest seq
    %in between, add 2n and 2n+1
    xp(slice,n/2-1) = abs(x(1));
    xp(slice,n/2-n/2+1) = abs(x(n));
    xp(slice,1+n/2-(2:n/2)) = abs(x(2:2:n-1))+abs(x(3:2:n-1));
    slice = slice + 1; 
end

subplot(2,2,1)
specgram(y,n,fs,[],0)
colormap(gray)
im = getframe;
im = sum(im.cdata,3) ;

subplot(2,2,2)
for i=1:size(im,2)
    [mm,mp]=max(im(:,i));
    imm(:,i) = 0;
    imm(mp,i) = 1;       
end
imagesc(imm)
%imagesc(im>0.5*max(max(im)))
colormap(gray)

subplot(2,2,3)
xp = (xp');
imagesc(log(xp))
colormap(gray)

subplot(2,2,4)
for i=1:size(xp,2)
    [mm,mp]=max(xp(:,i));
    xp(:,i) = 0;
    xp(mp,i) = 1;
end
imagesc((xp))

%imagesc(xp>0.5*max(max(xp)))
colormap(gray)