function [slope range] = evdnaslopec(dna, deg, step)
% takes a DNA sequence (length must be power of 2) and produces
% evolutionary slopes

N = length(dna);
coorda=(dna == 'A');
coordc=(dna == 'C');
coordg=(dna == 'G');
coordt=(dna == 'T');
null = zeros(size(coorda)); 
Y=[coorda; coordc; coordg; coordt; ...   
   coorda; coordg; coordc; coordt; ...
   coorda; coordt; coordc; coordg; ...
   null; null; null; null];
Y = cumsum(Y,2);

         %null's extend dimension to power of 2
haar=[sqrt(2)/2 sqrt(2)/2];
% haar1=[  -0.075765714789341  -0.029635527645954   0.497618667632458 ...
%    0.803738751805216   0.297857795605542  -0.099219543576935...
%   -0.012603967262261   0.032223100604071];

lener1=[]; lener2=[];
window = 2^deg;
stepsize = 2^step;
range = (window/2+1):stepsize:(N-window/2);
W2 = WavMat(haar,16);
Wk =  WavMat(haar,window);
D = W2 * Y;

for k = range
Zk = D(1:16,(k-window/2):(k+window/2-1)) * Wk';
lener1 =[lener1 log2(mean(mean(Zk(5:7, window/4+1:window/2).^2)))];
lener2 =[lener2 log2(mean(mean(Zk(9:14, window/2+1:window).^2)))];
end

slope = lener2 - lener1;
nn = length(slope);
close all
figure(2)
plot(range, slope,'k-');
hold on
xlim([0 max(range)])
%plot(range, -2 * ones(size(range)),'r-','linewidth',2)
%plot(range(1:nn/2), mean(slope(1:nn/2)) .* ones(size(range(1:nn/2))),'g-','linewidth',2)
%plot(range(nn/2:end), mean(slope(nn/2:end)) .* ones(size(range(nn/2:end))),'g-','linewidth',2)
end
